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CHAPTER  I 


INTRODUCTION 


Introdactioii 

Composite  materials  have  become  widely  used  in  engineering  applications  in 
the  past  coiq)le  of  decades.  This  class  of  materials  holds  many  benefits  ^^iien  used 
apfvopriately  in  engineering  applications.  Because  of  analysis  uncertainties  many 
composite  components  are  "over-engineered"  and  the  design  is  often  governed  by 
reiterative  component  testing.  In  tiiese  cases,  the  full  benefit  of  composite  materials  is 
not  realized.  This  has  led  to  the  develojHnent  of  analysis  aids  for  several  different 
structural  member  types. 

One  of  the  major  composite  structural  members  is  the  composite  plate.  A 
plate  is  a  load  carrying  member  ^ch  is  bounded  by  two  parallel  planes  called  feces. 
Each  face  has  the  same  characteristic  length  and  width  dimensions  and  are  bounded 
by  the  plate  edges.  The  distance  between  these  feces  is  the  plate  thickness  and  this 
thickness  is  considered  tc  «  sirudl  compared  to  the  dimensions  of  the  faces.  The 
plate  feces  can  take  on  many  different  types  of  shapes  (rectangular,  circular,  elliptical 
and  others).  Composite  plates  have  been  used  in  aeronautical  structures  for  years. 
Composite  plates  are  currently  being  used  in  kmd-based  construction  because  of  their 
exceptional  environmental  properties. 
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There  are  currently  analysis  im)grams  which  include  composite  materials.  In 
addition  to  specific  composite  {vograms,  several  general  finite-element  iMX>gr8ms 
available  on  large  systems  incorporate  the  analysis  of  composite  materials.  Analysis 
techniques  for  composites  have  been  changing  rapidly.  Since  these  larger  {xograms 
have  included  composites  as  an  auxiliaiy  component,  they  do  not  always  keep  up  with 
current  research  in  this  area.  Also,  the  size  of  these  programs  inohibit  tl^ir  use  on 
microcomputers.  Some  authors  have  publisl^  computer  jvograms  for  specific 
composite  structures  or  limited  composite  material  lay-iq)6.  No  author,  however,  has 
published  a  computer  program  for  the  analysis  of  general  composite  plates. 

Objectives 

The  main  objective  of  this  research  is  to  jncduce  a  woiking  computer  program 
for  the  analysis  of  general  composite  plates  to  be  used  on  microcomputers.  The 
{Hogram  jnesented  in  this  paper  is  limited  to  the  analysis  of  laminated  composite 
plates  with  elastic  behavior  and  small  deflections.  Shear  deformation  is  included  in 
the  analysis  because  of  the  material  behavior  response  specific  to  composites.  This 
program  is  an  revision  of  an  existing  jxcgram  published  by  J.  N.  Reddy  [25].  Reddy’s 
program  was  developed  to  analyze  orthotropic  materials  with  elastic  behavior  and 
snudl  deflections.  Although  single-layered  composites  exhibit  this  behavior,  most  of 
the  composite  plates  used  in  applications  have  more  than  one  layer  and  require  a  more 
complex  jnogram  for  analysis. 

To  validate  the  computer  code,  results  from  the  program  are  compared  against 
results  from  other  analytical  methods  and  results  from  other  authors  in  the  literature. 
This  test  is  used  to  insure  that  the  program  im)perly  employs  first-order  shear 
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deformation  theory  (FSDT)  of  composite  plates.  Additionally,  the  effect  of  shear 
deformation  in  composite  plates  is  observed  by  comparing  the  inx>gram  results  against 
results  from  a  plate  theory  that  does  not  include  shear  deformatioiL  Because 
comprehensive  instructions  and  documented  source  code  are  included,  the  program 
should  prove  to  be  a  valuable  educational  aid  for  teaching  the  application  of  the  finite 
element  method  to  composite  structures  in  advanced  composite  classes. 

Overview 

Chapter  n  begins  with  an  introduction  to  composite  materials  with  a 
background  on  their  mechanics.  An  introduction  to  current  composite  plate  theories 
ends  the  chapter.  Definitions  of  variables  and  sign  conventions  used  in  the  program 
and  the  rest  of  the  paper  are  presented  to  aid  in  the  reader's  comprehension.  Chapter 
in  provides  an  introduction  and  derivation  of  first-order  shear  deformation  theory  of 
composite  plates  using  variational  energy  formulation.  Chapter  IV  shows  how  this 
theory  is  transformed  into  a  finite  element  model  for  use  in  the  computer  program. 
Numerical  results  from  the  computer  program  are  compared  against  other  analytical 
method  solutions  to  validate  the  program  code  in  Chapter  V.  Finally,  conclusions 
derived  from  the  results  and  recommendations  for  future  woric  are  presented  in 
Chapter  VI. 

It  is  assumed  that  the  reader  has  a  general  knowledge  of  composite  materials, 
plate  theory,  and  the  finite  element  and  variational  methods.  Some  background  is 
presented  in  these  areas  to  define  terms  and  conventions  used  in  the  plate  theory.  For 
further  information  in  these  areas,  see  the  following  references:  composite  materials 
[1,30],  plate  theory  [27,29,33],  finite  element  and  variational  methods  [8,10,21,24,25]. 


CHAPTER  n 


MECHANICS  OF  COMPOSITE  MATERIALS  AND  PLATE  THEORIES 

Composite  Materials 

Th^  are  many  types  of  composite  materials  used  in  the  fabrication  of 

structural  components.  The  term  "composite"  refers  simply  to  a  material  made  of 

more  than  one  distinct  constituent  Composites  have  become  known  as  nuiterials 

Ariiich  have  clear  boundaries  between  its  constituents,  and  ^riiose  constituents  have 

markedly  different  material  properties.  The  constituents  combine  to  form  a  composite 

material  with  material  properties  considerably  different  from  any  of  its  constituents. 

Most  of  the  modem  composites  contain  either  particulates  or  fibers  as  main 

constituents.  Particle-reinforced  composites  are  formed  by  suspending  either  random 

or  preferred  orientation  particles  in  a  surrounding  material.  The  material  properties  of 

these  types  of  composites  are  obtained  finm  load  tests  and  are  similar  to  isotropic  (for 

random-oriented  particulates)  or  orthotropic  (for  preferred  orientation  particulates) 

materials.  Fiber-reinforced  composites  are  made  of  fibers  suspended  in  a  surrounding 

material.  The  fn>ers  may  be  either  continuous  or  discontinuous  (short-fiber).  Fiber- 

reinforced  composites  may  either  single-layered  (including  multiple  plies  of  the 

same  fiber  orientation)  or  multi-layered.  See  Figure  1  for  an  outline  of  composite 

classifications.  The  program  ixesented  in  this  paper,  COMPLATE,  is  useful  for 

analyzing  all  of  the  above  composite  types.  The  most  general  case  of  composites  are 

multi-layered  continuous-fiber-reinforced  hybrid  composites.  The  mechanical 
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response  of  the  other  composite  types  can  be  modeled  with  simplifications  to  this 
general  case.  Also,  continuous-fiber-reinforced  composites  are  the  most  commonly 
used  composites  for  structural  components  where  hi^  strength  is  required  For  these 
reasons,  the  mechanics  of  continuous-fiber-reinforced  composites  are  defined  further 
and  are  utilized  in  the  development  of  the  computer  program. 


Coiiq)osite  materials 


Fiber-retoforced  composites  Partide-reinfoTced  composites 

(fibrous  composites)  (particulate  composites) 


RandcMii  Prefierred 
orientation  orientation 

I - - 1 

Single-layer  con^wsites  Multilayered  composites 
(laminae)  (laminates) 


l.aminates  l^brida 

I  I 

Continuous-fiber-reinforced  Discontinuous-fiber-reinforced 
Gonq)osites  composites 

I - ' - 1  I - ' - 1 

Unidirectional  Bidirectional  Random  Preferred 

reinfbroement  reinforcement  orientation  orientation 

(wov«i) 

Figure  1;  Types  of  Composite  Materials  [1]. 

Continuous-fiber-reinforced  composites  (hereafter,  simply  composites)  differ 
in  many  ways  from  isotropic  materials.  Composites  are  generally  composed  of  two 
distinct  materials:  reinforcements  (fibers)  and  matrix  (bonding  material). 
Reinforcements  made  of  fibers  form  the  strength  of  a  composite  because  they  cany  a 
majority  of  the  load  Matrix  is  the  material  in  between  these  fibers  that  binds  the 
fibers  and  provides  for  load  transfer  between  fibers  in  case  of  fiber  breakage.  The 
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matrix  also  {votects  the  fibers  from  environmental  degradation  and  damage  due  to 
handling.  The  matrix  nuderial  generally  has  strength  and  stiffoess  jxoperties  much 
less  than  the  reinforcement  or  fiber. 

Composites  often  achieve  stiength-to-weight  ratios  significantly  higher  than 
metals.  Atomic  theory  predicts  strengths  much  higher  than  those  actually  found  in 
I»actice  for  all  materials.  The  reason  for  this  shortfall  in  strength  arises  from  inherent 
defects  at  both  the  microscopic  (atomic)  and  macroscopic  (visible)  levels  created 
during  material  processing.  The  lai^est  allowable  defect  size  at  the  nuu^roscopic  level 
depends  on  the  cross-sectional  area  of  the  material.  For  bulk  materials,  relatively 
large  defects  can  occur  during  material  inocessing.  For  the  noanufacture  of  composite 
fibers,  the  size  of  defects  is  reduced  because  the  cross-sectional  area  of  the  fiber  is 
relatively  small.  If  a  visible  defect  is  fxesent  in  the  fiber  material,  it  breaks  as  it  is 
stretched  during  manufacture.  The  unbroken  portion  of  fibers  have  defect  sizes 
limited  to  the  microscopic  level.  By  themselves,  fibers  are  not  useful  for  structural 
applications  because  of  their  small  size  and  strength.  A  large  number  of  fibers  are 
bonded  together  by  use  of  a  matrix  to  form  a  hi^stiength  material.  There  are  many 
methods  for  manufacturing  composite  materials.  The  main  concern  of  this  paper  is 
composite  plate  q>plications,  so  tl^  following  discussion  refers  to  the  structure  of 
composite  plates.  However,  for  more  information  pertaining  to  the  manufacturing  of 
composites  see  references  by  Agarwal  and  Broutman  [1],  and  Vinson  and  Sierakowski 
[30]. 
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Lamina  and  Laminata 

The  ply  is  tl^  basic  building  block  of  composite  plates.  A  ply  is  the  thin  sheet 
of  unidirectional  fibers  bonded  by  matrix  material  developed  during  manufacture. 
The  ply  is  often  many  fiber  diameters  thicL  A  lamina  or  layer  is  formed  wdren  a 
unidirectional  ply  or  combination  of  unidirectional  plies  of  the  same  material  with  the 
same  global  fiber  orientation  is  suspended  in  a  matrix.  Although  it  may  consist  of 
several  plies,  the  important  aspect  of  the  lamina  is  that  it  is  defined  as  a  layer  of 
material  with  common  directional  material  {xoperties.  A  multi-plied  lamina  contains 
a  fiberless  interface  between  plies  ^^ch  is  relatively  thin  and  is  often  ignored  for 
analysis  purposes.  In  practice,  fibers  are  not  equally  spaced,  but  for  schematic 
purposes,  the  lamina  is  often  depicted  having  a  single  layer  of  fibers  with  universal 
fiber  spacing  as  in  Figure  2. 

The  material  (x^operties  of  composites  differ  from  isotropic  materials  in  the 
following  way.  Each  lamina  exhibits  a  generalized  orthotropic  behavior  whose 
properties  are  different  on  three  mutually  perpendicular  planes  aligned  with  the  fiber 
direction  shown  as  1,  2,  3  in  Figure  2.  Material  properties  are  defined  in  the  three 
directions  corresponding  to  these  planes. 


Figure  2:  Schematic  view  of  a  lamina 
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These  lamina  are  then  stacked  with  the  fibers  aligned  at  different  angles  to 
form  what  is  called  a  laminate  as  shown  in  Figure  3.  The  lamina  are  labeled 
according  to  their  fiber  angle  relative  to  a  global  direction  (x-axis).  A  code  has  been 
developed  to  label  laminate  stacking  sequences.  For  example,  [0/45/90]  is  a  laminate 
composed  of  three  equally  thick  lamina  whose  fibers  are  oriented  0**,  45**,  and  90° 
respectively  to  the  principle  reference  direction  starting  with  the  bottom  layer  (as 
shown  in  Figure  3).  A  subscript  s,  [90/45/0]^,  denotes  a  symmetric  lay-iq)  ^^liere  the 
top  layers  are  stacked  in  reverse  order  or  [90/45/0/0/45/90]  and  a  numerical  subscript 
denotes  the  number  of  repeated  plies,  [902/45y0J  =  [90/90/45/45/45/45/0/0]  for 
example. 

As  the  lamina  are  stacked  to  form  a  laminate,  effective  macroscopic  ivoperdes 
are  developed  to  characterize  the  laminate.  These  inoperties  are  assumed  to  be 
homogeneous  although  direction  dependent  (anisotropic)  and  are  a  weighted  average 
of  the  properties  of  the  composite  constituents.  Therefore,  two  laminates  made  of  the 
same  fiber  and  matrix  ’^^^terial  may  have  very  different  macroscopic  material 
fvoperties  because  of  a  difference  in  their  stacking  sequence. 


Figure  3:  Schematic  of  a  three-layered  laminate  [0/45/90]. 
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General  strengths  and  properties  are  experimentally  determined  for 
unidirectional  lamina  of  composite  materials.  These  lamina  strengths  are 
incorporated  through  the  use  of  equations  to  predict  effective  macroscopic  {voperties. 
Strengths  of  laminate  stacking  sequences  are  determined  by  one  of  several  failure 
theories  [1].  Test  specimen  are  used  to  experimentally  measure  material  properties 
and  consist  of  small  strips  of  composite.  These  specimen  are  checked  for  af^nrent 
flaws  or  defects  and  their  edges  are  smoothed.  The  test  specimen,  dierefore,  form  an 
ideal  base-line  on  the  strength  of  the  composite. 

Lamina  Constitutive  Relations 

For  a  given  lamina,  the  stiffness  properties  are  generally  given  with  respect  to 
inincipal  fiber  directions.  Direction-1  is  aligned  with  the  longitudinal  direction  of  the 
fibers.  Direction-2  is  aligned  with  the  direction  transverse  to  the  fibers  in  the  lamina 
plane.  Direction-3  is  normal  to  both  the  1  and  2  directions.  Figure  4  shows  these 
directions  with  respect  to  fiber  alignment 
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1 


denotes  lamina  fiber  direction 
Figure  4:  Principal  fiber  directions  (1,2,3). 

The  material  properties  are  determined  throu^  experimentation.  Most  lamina 
a^e  characterized  by  the  following  independent  material  properties:  E„  Ej,  v,2,  Gjj, 
G23,  G,3  .  These  properties  are  used  to  develop  the  stiffiiess  matrix. 

A  general  6x6  orthotropic  stiffness  matrix  relates  the  6  principal  normal  and 
shear  strains  to  the  corresponding  principal  stresses  [1].  For  the  laminate  plate  theory 
presented  in  this  paper,  the  out-of-plane  normal  strain,  e, ,  is  assumed  to  be  zero.  This 
strain  is  uncoupled  from  the  other  strains  and  it  allows  the  stifrhess  matrix  to  be 
reduced  to  5  x  5.  For  each  lamina,  the  orthotropic  stiffness  nutirix  aligned  with 
principal  fiber  directions  is  defined  by  the  following  stress-strain  relationship  given  in 
equations  2-1  and  2-2. 


Qn  Qi2  0  0  0  C] 

Qi2  Q22  0  0  0  62 

0  0  Qj3  0  0  <  y,2  • 

0  0  0  Q44  0  723 

_  0  0  0  0  Q55.  .Yi3. 


(2-1) 


>^4iere  the  matrix  terms  are  defined  as: 
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Qn=  — 


l-V,2Vj, 


l-VuVj, 


Q33  =  Gi2 

J  ^12^2 

and  V,,  =  — 

E. 


Q44  ~  ^23 


Q22  = 

Q55  ^13 


^  ''12 ''21 


(2-2) 


The  stiffiiess  matrix  given  above  is  most  useful  for  characterizing  lamina 
properties.  Laminates  are  formed  by  stacking  layers  of  lamina  with  varying  fHrer 
orientations,  thicknesses,  and  materials.  To  accommodate  the  variance  in  fiber 
orientation,  the  lamina  stiffness  matrices  must  be  transformed  to  a  common  global 
orientation.  For  each  lamina  the  fiber  orientation  is  defined  by  the  angle,  6^ ,  that  the 
fiber  direction  makes  with  the  xnaxis.  The  angle  6,;  is  defined  as  positive  in  the 
counterclockwise  direction  and  negative  in  clockwise  direction  as  shown  in  Figure  5. 
The  angle  0,^  can  have  a  value  between  90**  and  -90**. 


Figures:  Global  fiber  orientation,  0. 
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The  following  relationships  are  based  on  a  laminate  of  N  layers.  For  each 
lamina  (k  =  1^,3,...>0>  transformed  stiffoess  matrix,  [Q],  is  defined  the  stress- 
strain  relationship  aligned  with  the  global  axes  (x,y,z)  as  follows: 


(k) 

Qn 

Q,2 

Q.3 

0 

0  ■ 

(k) 

•  ' 

e. 

Q.2 

Q22 

Q23 

0 

0 

>  = 

Qi3 

Q23 

Q33 

0 

0 

i 

y^y 

0 

0 

0 

Q44 

Q45 

Yyx 

0 

0 

0 

Q45 

Q53. 

y  X* . 

(2-3) 


where  for  each  lamina,  k: 

Qii  “  Qii™  ■1'2(Q,2 +2Q33)m  n^+Qjjn^ 

^12  ~  (Qn  ■^Q22  “  ®  ) 

Qi3  =  Qiim’n-Q22™®’  -(Qi2  +2Q33)mn(m^  -n*) 

Q22  ~  Qu®  2(Qi2  2Q33)m  n*  +Q22na 

^23  =  Q„mn’  -Q22m’n  +  (Q,2  +2Q33)mn(m*  -  n*)  (2-4) 

^33  =  (Qii  +Q22  -  2Q,2  )m*n*  +Q33(m*  -  n*  f 
Q-m  =  Q44®*  +Q3J®* 

Q4s=(Q35-Q«)“n 

Qjj  =  Q^^n^ +Q,jm* 

and  m  =  cosO^,  n  =  sinO^. 
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Lamiiiate  Constitative  Relations 

In  order  to  develop  constitutive  relationships  that  are  independent  of  z,  it  is 
useful  to  define  load  and  moment  resultants.  These  resultants  are  the  loads  and 
moments  per  unit  length  along  the  lamina  x  and  y  cross-sections  and  acting  through 
the  laminate  mid-plane.  The  orientation  and  positive  direction  are  depicted  in  Figure 
6. 


The  in-plane  load  resultants  (N, ,  Ny ,  )  are  defined  as  the  integrals  of  the  in¬ 

plane  stresses  throu^  the  thickness  in  the  respective  directions  shown  above  in  Figure 
6. 
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N.=jlo.di  N,  =  jio,dz  N.,=j|t.,dz  (2-5) 

2  2  2 

The  moment  resultants  (M, ,  My .  M^)  are  defined  as  the  integrals  of  tte 
moments  created  about  the  laminate  mid-plane  by  die  in-plane  stresses  through  the 
thickness  in  the  respective  directions  shown  above  in  Figure  6. 

M,  =  Jio,zdz  My  =  JiOyZdz  M^  =  j4T^zdz  (2-6) 

2  2  2 

The  shear  load  resultants  (Q,  and  Qy)  are  defined  as  the  integrals  of  the 
transverse  shear  stresses  throu^  the  laminate  thickness. 

Qx  =  JA  Qy  = 

2  2 

Figure  7  shows  the  laminate  coordinate  system  throu^  the  laminate  thickness 
with  terms  used  in  the  following  sets  of  equations. 


z 


15 


Figure  7;  T  iiminate  coordinates  and  terms. 


The  load  resultant  equations  2-5  to  2-7  along  with  the  stress-strain  relations  2- 
3  and  2-4  are  used  to  develop  the  laminate  stiffness  matrices.  The  extensional 
stiffness,  coupling  stiffness,  and  bending  stiffness  matrices,  A^,  D^,  respectively, 

are  defined  by  the  following  matrix  equations: 


(2-8) 


(2-9) 


)6 


where  are  the  in-plane  strains  at  the  laminate  mid-plane,  K,,Ky,K^  are  the 

laminate  curvatures,  and  are  the  average  transverse  shear  strains  as  defined 

in  the  next  chapter  by  first-order  shear  (^formation  theory. 

The  defmitions  of  A,  B,  D  plate  stiffness  matrix  terms  (i,  j  =  1,2,3)  are  given 
followed  by  simplifications  for  laminated  plates. 


A,  =p,QS'>dz  = 

2  k-l  k-l 


B,  =  Ji(5S“ab  =  llQf  -  zi.,) 

k-I 


(2-10) 


2 

h 


k-l 


D, = jiQi‘>z‘dz = it5r(z:  -  zi., ) 

2  J  k-l  k-l 


/  .3\ 

tkZk  +  — 

'  ‘  12  y 


The  definitions  of  the  terms  z,  z^ ,  z^_^  .t^  and  z^  are  given  in  Figure  7.  The  transverse 
shear  matrix  terms  (i,  j  =  4,  S)  are  given  by  the  following  equations. 

^  N  N 

A,  =  k.JiQ?>dz=  k.2;Q®[z.  -z._,]=  (2-») 

2  k-l  k-l 

Mdiere  1^  is  the  shear  correction  factor.  Methods  for  determining  the  value  of  the 
shear  correction  factor  are  {vesented  in  [5,31,33,34] 

Inertial  properties  of  the  composite  plate  are  required  for  dynamic  analysis. 
The  density  of  each  lamina  is  given  by  p<^>  where  k  is  the  layer  number.  The  following 
terms  represent  transverse,  transverse-rotation  coupling,  and  rotational  inertial 
properties  respectively. 
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Pi  =  -  Zk_,  )  =  2p‘^Hk 

2  k-1  k-l 


P2  =  Jjp'‘^2dz  =  i5]p<^>(Z^-Z*_,)  = 

2  ^  k-l  k-l 


(2-12) 


P3  =  Jjp^^^z^dz  =  -  A-i )  =  Zp 

2  ■J  k-l 


N 


(k) 


k-l 


tk^k  +  — 

^  ^  12, 


The  relationships  for  general  orthotropic-layered  laminated  plates  given  above 
can  be  applied  to  many  other  types  of  composite  plates.  Single-layered  laminates  can 
be  accommodated  by  using  a  single  layer  (N=l)  in  the  above  equations.  For  preferred 
orientation  (both  fKirous  and  particulate)  and  bi-directional  lamina,  the  orthotropic 
relationships  along  with  the  preferred  orientation  (9)  can  be  ^lied  to  a  single  layer. 
For  random  orientation  lamina  (discontinuous-fiber  and  particulate),  isotropic 
relationships  are  obtained  by  using  a  single  modulus  of  elasticity  (E,  =  Ej  =  E)  a  single 
shear  modulus  (G,2  =033=  0,3  =  G),  and  Poisson's  ratio  (v,2  =  v)  in  the  above 
relationships. 

Plate  Theories 

Plate  theories  are  simplifications  to  general  three-dimensional  elasticity  theory 
and  were  developed  to  analyze  one  of  the  basic  structural  member  types,  the  plate. 
Three-dimensional  elasticity  theory  may  theoretically  be  used  to  analyze  any  solid 
object  Elasticity  theory  is  often  prohibitive  in  practical  use  because  of  the  complexity 
of  the  solid  object  and  the  cost  of  aj^lying  the  general  theory  to  each  analysis  case. 
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Because  of  the  general  shape  and  load  conditions  common  to  plates,  certain 
assumptions  are  made  which  reduce  the  equations  of  elasticity  theory  to  several  less- 
complex  governing  equations.  There  are  two  important  points  to  remember  vdien 
iqrplying  plate  theories.  First,  the  assumptions  mark  in  simplifying  the  ^veming 
equations  limit  the  types  of  cases  ^^re  a  particular  plate  theory  can  be  effectively 
ai^lied.  Second,  the  assumptions  create  an  approximate  solution  for  the  plate 
problem.  Other  methods  including  testing  should  be  used  to  determine  the  accuracy 
of  the  plate  theory  solutiorL  Many  of  the  plate  theories  were  first  developed  for 
isotropic  materials  and  were  later  adapted  to  composite  materials. 

Two-dimensional  plate  theories  are  developed  by  assuming  displacement 
fimctions.  These  fimctions  are  characterized  by  equation  2-13.  The  displacement  of 
any  point  in  the  plate  (u,  ,U2 ,03)  is  defined  by  its  mid-plane  displacement  (u,v,w),  a 
fimction  of  the  mid-plane  coordinate  (x,y),  and  an  assumed  form  of  the  displacement 
through  the  plate  thickness  (U,  V,W),  dependent  on  the  laminate  mid-plane  coordinate 
(x,y)  and  distance  from  the  mid-plane  (z).  These  fimctions  are  also  dependent  on  time 
(t)  for  the  case  of  dynamic  problems.  By  separating  the  displacement  functions  into 
these  two  parts,  the  analysis  can  be  reduced  from  three  dimensions  (x,y,z)  to  two 
dimensions  (x,y). 

Ui  (x,y,z,t)  =  u(x,y,t)  +  U(x,y,z,t) 

U2(x,y,z,t)  =  v(x,y,t)  + V(x,y,z,t)  (2-13) 

U3(x,y,z,t)  =  w(x,y,t)+W(x,y,z,t) 

Classical  plate  theory  (CPT)  was  the  first  form  of  the  plate  theories  and  is 
attributed  to  Kirkhhoff  [20].  This  theory  was  the  first  attempt  to  characterize  thin 
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isotropic  plates  and  is  limited  in  its  scope  of  plications  due  to  many  assumptions. 
The  adaptation  of  this  theory  to  laminated  composite  materials  is  generally  attributed 
to  Yang,  Norris  and  Stavsky  with  additions  by  Whitney  and  Pagano  [33].  The 
following  assumptions  provide  the  basis  for  CPT  [35]. 

1 .  Plane  sections  of  the  plate  cross  section  remain  plane  and  normal  to  the 
mid-surface. 

2.  The  deflections  are  small  compared  to  the  plate  thickness. 

3.  Transverse  normal  strain  is  zero  and  transverse  shear  strains  are  negligible. 

4.  Transverse  normal  stress  is  negligible. 


These  assumptions  result  in  the  simplified  displacement  functions  given  in 
equation  2-14.  These  functions  have  three  degrees  of  fieedom  (u,v,w)  which  are 
dependent  on  x,  y,  t  only. 


dw 

u,  (x,y,z,t)  =  u(x,y,t)  -  z—(x,y,t) 

dx. 

dw 

U2(x,y,z,t)  =  v(x,y,t)-  z~(x,y,t) 

dy 


(2-14) 


U3(x,y,z,t)  =  w(x,y,t) 


This  theory  is  adequate  for  a  large  class  of  isotropic  plates  and  very  thin 
composite  plates.  For  thicker  plates,  with  length  and  width  to  thickness  ratios  of  less 
than  10,  CPT  tends  to  under  jxedict  the  deflections  in  plate.  This  is  caused  from  the 
transverse  shear  strains  being  larger  than  the  assumption  requires.  This  i»oblem  is 
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also  encountered  in  moderately  thin  composite  plates  because  of  the  directional 
material  jvoperties  unique  to  composites. 

One  theory  that  addresses  this  inadequacy  is  the  first-order  shear  deformation 
theoiy  (FSDT).  This  theory  was  developed  by  Reissner  for  static  analysis  and  refined 
by  Mindlin  for  (fynamic  analysis  of  isotropic  plates  [18].  Yang  et  al.  modified  this 
theoiy  for  composite  plates  with  further  refinements  by  Whitney  and  Pagano  [33]. 
This  theoiy  includes  transverse  shear  strain  in  the  analysis  and  gives  better  results  for 
deflections  and  stresses  in  composite  plates.  Equation  2-1 S  shows  the  assumed 

displacement  fimctions  for  this  theoiy.  Note  that  this  theory  allows  five  degrees  of 
freedom  (u,v,w,\t>,,\|;y).  This  tteoiy  is  known  as  first-order  because  the  total 

displacements  are  assumed  to  be  linear  fimctions  of  z  throu^  the  plate  thickness. 
FSDT  is  presented  in  greater  detail  in  the  next  chiqjter. 

u,  (x,y,z,t)  =  u(x,y,t)  +  zv,(x,y,t) 

U2(x,y,z,t)  =  v(x,y,t)  +  zVy(x,y,t)  (2-15) 

U3(x,y,z,t)  =  w(x,y,t) 


FSDT  does  not  adequately  address  boundary  conditions  on  the  plate  faces  or 
predict  the  interlaminar  shear  stresses  through  the  plate  thickness.  In  response, 
several  higher-order  shear  deformation  theories  (HSDTs)  have  been  fxesented  based 
on  work  by  Reissner  and  Schmidt  for  isotropic  plates  [18].  HSDT  has  been  extended 
to  laminated  plates  by  Nelson  and  Lorch,  Librescu,  Lo  et  al.  and  Red^  [20],  These 
hi^er-order  theories  retain  higher-order  terms  of  z  in  the  displacement  flmction  (see 
equation  3-1)  than  in  CPT  or  FSDT.  Several  HSDTs  have  been  jnesented  by  a 
number  of  authors.  A  sample  of  these  theories  can  be  found  the  following  references 
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[3,6,12,13,14,18,19,22,26^8,35].  Equation  2-16  shows  the  general  form  of  HSDPs. 
The  number  of  degrees  of  freedom  depends  on  the  order  of  the  displacement  functions 
in  terms  of  z  and  the  assumptions  of  the  particular  higher-order  theory. 

u,  (x,y,2,t)  =  u(x.y,t)  +  zv,  (x,y.t)+ z*C  Jx,y,t)  +  z’*,  (x,y,t)+-  • 

Uj  (x,y,z,t)  =  v(x,y  ,t)  +  z\(iy  (x,y,t)  +  z*Cy  (x,y,t)  +  z’^,  (x,y,  t)+-  •  •  (2-16) 

U3(x,y,z,t)  =  w(x,y.t)  +  ZM/,(x,y.t)  +  z*C*(x.y.t)+-- 


MSDrs  produce  more  accurate  transverse  shear  stress  results  than  the  two 
previous  theories,  however  the  resulting  deflection  and  normal  stresses  show  little 
improvement  over  FSDT  [22].  They  also  require  large  mathematical  and 
programming  costs.  Developing  a  program  for  one  of  these  theories  for  a 
microcomputer  is  prohibited  by  the  current  computing  ciqjacity  of  existing 
nucrocomputers.  FSDT  is  chosen  for  use  in  the  computer  program  for  its  relatively 
small  computing  requirements  balanced  with  its  improved  analysis  results. 


CHAPTERm 


FIRST-ORDER  SHEAR  DEFORMATION 
THEORY  OF  COMPOSITE  PLATES 


General  theory 

The  ivst-order  shear  deformation  theory  (FSDT)  presented  in  this  chapter  is 
used  for  the  revision  of  the  computer  i»Dgram  (COMPLATE)  from  the  {deviously 
written  fxogram  PLATE  by  J.  N.  Redcfy  [25].  This  theory  models  general  laminated 
composite  plates  and  includes  (fynamic  considerations.  The  main  purpose  of  utilizing 
this  theory  is  to  transform  a  three-dimensional  elasticity  problem  into  a  two- 
dimensional  problem.  The  energy  formulation  is  used  to  generate  mass  and  stiffness 
matrices  for  application  in  the  finite  element  method  jaesented  in  the  next  chiqrter. 

The  need  for  this  theory  arises  from  the  invalidity  of  neglecting  transverse 
shear  deformation  in  CPT.  Transverse  shear  is  ix>  longer  negligible  in  diick  (dates  of 
length  to  thickness  ratios  less  than  10  for  isotropic  plates.  Also,  shear  defommtion  is 
significant  in  composites  with  length  to  thickness  ratios  much  larger  than  10.  This  is 
due  to  the  effective  elastic  modulus  along  the  fiber  direction  (E,)  being  much  larger 
than  the  transverse  shear  moduli  (G,,  ,  G23),  sometimes  by  the  order  of  25  to  40 
compared  to  2.6  for  a  representative  isotropic  material  [22]. 

The  assumptions  listed  in  the  next  section  describe  tbe  restrictions  to  this 
model  and  should  be  considered  ^en  using  this  program  for  engineering 
a{^lications.  The  global  axes  described  in  the  next  section  are  shown  in  Figure  8  with 
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an  example  rectangular  plate.  Note  that  the  xy-plane  corresponds  to  the  laminate 
midplane  and  the  z-coordinate  describes  tl^  distance  and  direction  (upward  or 
downward)  of  a  point  with  respect  to  this  plane. 


Figure  8:  Schematic  of  a  rectangular  plate. 


Assumptions 


The  basic  assumptions  for  this  first-order  shear  deformation  theory  (FSDT)  are 
given  as  follows.  The  terms  used  in  the  assumptions  are  described  on  the  following 
pages  [33]. 

1.  The  plate  is  constructed  of  an  arbitrary  number  of  orthotropic  layers 
(laminae)  which  are  perfectly  boned  together.  The  directions  of  principle 
orthotropic  material  symmetry,  the  thickness,  and  the  material  of  each  layer 
may  vary. 

2.  The  plate  is  considered  to  be  relatively  thin  compared  to  face  dimensions. 

3.  Plate  displacements  (u,  v,  w)  are  small  compared  to  the  plate  thickness  (h). 

4.  In-plane  strains  (e,,  Gy,  Y^)  are  small. 
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5.  To  include  in-plane  force  effects,  non-linear  terms  in  the  equations  of 
motion  involving  products  of  stresses  and  plate  slopes  are  retained  and  all 
other  non-linear  terms  are  neglected 

6.  Transverse  shear  strain  Yh)  are  included  in  the  analysis  in  case  they  are 
not  negligible. 

7.  The  total  in-plane  displacements  (u„  Uj)  are  linear  functions  of  the  z- 
coordinate  through  the  plate  thickness. 

8.  Transverse  normal  strain  (e,)  is  i^gligible  compared  to  other  strains. 

9.  Each  lamina  behaves  in  an  elastic  manner  and  is  governed  by  Hooke's  law. 

10.  The  total  plate  thickness  is  uniform  throughout  the  plate. 

11.  Body  forces  are  negligible  compared  to  other  plate  forces. 

12.  All  linear  inertial  terms  are  retained  for  dynamic  analysis. 


Variational  Energy  Formulation 

For  FSDT,  a  first-order  (linear)  displacement  Held  in  terms  of  z  is  assumed 
The  general  displacement  of  any  point  in  the  plate  is  described  by  the  following  first- 
order  displacement  fimctions. 

u,  (x,y,  z,  t)  =  u(x,y,  t)  +  z\|/,  (x,y,  t) 

U2(x,y,z,t)  =  v(x,y,t)  +  zVy(x,y,t)  (3-1) 

U3(x,y,z,t)  =  w(x,y,t) 


where  u, ,  Uj,  and  u,  are  the  displacements  in  the  x,  y,  and  z  directions  respectively,  u, 
V,  and  w  are  the  displacements  of  tlm  laminate  mid-plane  in  the  same  directions,  and 
Vx  and  %  are  the  rotations  in  the  xz  and  yz  planes  respectively  caused  by  plate 
bending  and  transverse  shear  defoimatioa  By  adhering  to  the  assumption  of  small 
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displacements  the  strains  are  derived  from  the  displacement  functions  in  equation  3-1 
as  follows; 


The  definitions  of  the  mid-plane  strains  (e°,e°,Y^)  and  the  curvatures 
(K,,Ky,K^)  can  easily  be  derived  from  die  last  two  equalities  in  each  equation. 
Note  that  the  transverse  shear  strains  (Y^,Yn)  are  constant  through  the  plate 

thickness  (independent  of  z).  Since  achial  transverse  shear  strains  are  not  constant 
through  the  plate  thickness,  the  ones  {nedicted  by  this  theory  refHesent  average  shear 
strains. 

This  FSDT  is  based  on  a  displacement  derivation.  For  each  point  in  the 
laminate  mid-idane,  five  degrees  of  freedom  are  defined.  This  allows  for  enou^ 
degrees  of  freedom  to  provide  adequate  results  for  a  majority  of  composite  plate 
applications.  These  degrees  of  freedom  are  also  referred  to  as  generalized 
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displacements  and  consist  of  tlM  following  disfdacements  and  rotations;  u,  v,  w,  y, 
and  .  Each  generalized  displacement  is  a  function  only  of  mid-plane  position  (x,  y) 
and  time  (t)  for  dynamic  cases.  These  generalized  displacements  allow  for  the 
reduction  of  three-dimensional  model  (x,y^)  to  a  two-dimensional  mcxiel  (x,y)  of  the 
laminate  mid-plane  for  analysis  purposes. 

For  the  general  dynamic  case,  tlM  energy  formulation  is  based  on  Hamilton's 
Principle: 

f8Ldt  =  0  (3-3) 

The  Lagrangian  (L)  may  be  defined  as  components  of  energy  in  the  following 
manner. 


L  =  U  +  V-T  (3-4) 

where  U  is  the  total  strain  energy,  V  is  the  potential  energy  due  to  the  uniformly 
distributed  transverse  load,  and  T  is  the  total  kinetic  energy  of  the  plate. 

The  fhst  variation  of  the  Lagrangian  can  be  written  as  the  variation  of  its 
components: 

5L  =  5U  +  5V-5T  (3-5) 

In  order  to  fmd  the  first  variation  of  the  Lagrangian,  the  first  variation  of  each 
of  the  components  are  derived  in  terms  of  displacements.  The  total  strain  energy  of 
the  plate  (U)  is  defined  by  the  integration  of  the  strain  energy  in  terms  of  stresses  and 
strains  of  each  point  in  the  plate. 


27 


u  =  iJJJ(oxex  +OyE,  +o,e,  +x.yY,y  +T^Y^  +T„Y„)dV  (3-6) 

V 


Using  the  definitions  of  stress  from  equation  2-3  in  the  fvevious  chapter,  the 
strain  energy  can  be  recast  as: 


U  -  i’ Q22®y  ■*'Q33Yxy  ‘•'^Qu^x^y  "*■  2Ql3®xYxy  ■*‘2Qi2Ex6, 

V 

+2Qj3eyYxy  +Q44Y^  +Q5jYL  +2Q«Y,xY„]dV 


Substituting  the  definitions  of  strain  from  equation  3-2,  the  following  equation 
is  obtained: 


Expanding  the  equation  and  factoring  terms  of  z  (1,  z,  z^)  yields  the  following: 
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(3-9) 


Since  u,  v,  w,  and  \\fy  are  independent  of  z,  the  triple  integral  may  be 
integrated  with  respect  to  z,  and  the  relations  for  and  from  equations  2-10 
and  2-1 1  may  be  iqiplied  to  yield; 


U=ir/A„^'+2B,.^^+D„i5£-^*  avSv,  av, 

"ax  "ax  ax  "  ^ 


'L 


+A 


33 


ax 


2  ^2 

-  +A2,—  +2B2, - 

dy  dy  dy 


+  D, 


dy 


au  ^  av 

dy  dx 


+  2B 


'"au  av^^ 


33 


dy  dx)\  dy  dx 


+  D 


33 


a»|/,  ^  a\|/. 


- .  au  av  _ 
■*■2^12  -  -  ■*■^12 


dx  dy 


au  ^  av  ay, 
dx  dy  dy  dx 


+  D  I  2A 

^^12  -  a.  ^^^13 


dx  dy 


dy  dx 

f 
\ 


y 

> 


au  au  ^  au  av 

dxdy  dxdx 


+2B 


13 


'au 

fav,  1  s^^fy^ 

du  dv^ 

av 

1  ‘'“X 

"aVx  ,  ^y' 

^dx 

^  dy  dx  j 

dx 

.ay'^ax^ 

^  ay  ax  ^ 

30 


n  5 


dy  dx.  )\dy 


5u  ^  „  I  5u  dw 


— (8m/,)+  B,3— +Bj3  — +  B33  — +  —  +D,3 


dx  dy  I8y  dx 


au/,  f 0«i,  a  a  f 

(  5w)l  5  ,  r  -  Z'  Sw')  f  ^  /-R  \l  JA 


(3-11) 


Using  the  expression  for  the  derivative  of  a  product,  the  following  relation  is 


derived: 


-l-m+f^±m]dA 

dx  dy  J 


'ft 

dy 

•JjfA(f,5u)  +  |- 

K^dx  dy 


(fjSu)  dA 


(3-12) 


Green's  Theorem  is  used  to  evaluate  the  last  integral  above: 


fff— (f,8u)  +  — (f28u)  dA=  f(f,8udy-f28udx) 

aV^  ay  ;  J 


(3-13) 


This  leads  to  the  following  equation  used  to  evaluate  8U: 


jj|f,  A(8u)+f,|-(Su)jdA  =  -Jj|^^+^j8udA  +  |(f,8ndy-fj5udx)  (3-14) 
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5U  may  be  divided  into  tenns  with  common  factors  of 
displacement,  5U^  to  SU^,  and  a  boundary  term,  SUg. 

8U  =  SUfc  +8Us,  +8Uj,  +80*,^  +8Us 


Thus, 
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(3-15) 


(3-16) 


(3-17) 
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8U^  =  -JJ 
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f  I  I  2  ' 
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(3-20) 


By  utilizing  Green's  Theorem  and  flying  the  laminate  load  resultant 
relations  from  equations  2-8  and  2-9,  the  boundary  condition  component  becomes: 
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6Us  =  J(n,6u  +  N,,8v+Q,8w  +  M,5i(/,  +M,y8v/y)dy 
Sx 

(3-21) 

-  J  (n^8u  +  Ny8v  +  Q,8w  +  M,y8n/,  +  My8M/y  )dx 

St 


For  the  computer  program,  the  loading  on  the  plate  is  assumed  to  be  a  uniform 
pressure,  thus  the  potential  energy  in  the  plate  due  to  applied  transverse  pressure  (V) 
is: 


V  =  JJqwdA  (3-22) 

A 

and  the  first  variation  of  V  is 

8V  =  JJq8wdA  (3-23) 

A 

The  total  kinetic  energy  of  the  plate  (T)  is  the  final  component  of  the 
Lagtangian.  It  consists  of  the  following  integral  of  the  energy  of  each  point  in  the 
plate. 


,  ^2^  ,  ^3' 
dt  dt 


dV 


(3-24) 


The  time  derivatives  of  the  displacements  are  found  by  taking  the  first  time 
derivative  of  equation  3-1  and  result  as  follows: 
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du,  du 
dt  a  dt 

du.  d\ 
dt  dt  dt 

aij  _  dw 

ir"ir 


(3-25) 


Substitutioii  of  these  time  derivatives  into  equation  3-24  yields: 


da^  dv^  aw"  -  f 
—  + —  + —  +2z 
dt  dt  dt 


gu  ^  dy 
dt  dt  dt  dt  ) 


+  z 


dt  dt 


2^ 


dV 


(3-26) 


Since  the  time  derivatives  of  u,  v,  w,  M't  and  are  independent  of  z,  the 
equation  may  be  integrated  with  respect  to  z  and  the  inertial  terms  from  equation  2-12 
(p, ,  p2 ,  and  P3)  may  be  applied: 


Pi 


—  +—  I 

dt  dt  dt 
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(3-27) 


Taking  the  first  variation  of  the  kinetic  energy  yields: 


— (8u)  +  p,-^— (8w) 


f  av  ^ 
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J^(8v.)+ 
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(3-28) 


dA 
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and  using  the  same  relation  as  the  one  derived  for  strain  energy  (equation  3-12),  5T 
becomes: 


5T-JJ 
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Pl^  +  P2-^ 


d^u  V 

'’=aF*'’>-aFr- 
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dA 


(3-29) 


The  above  variational  components  of  virtual  energy  (5U,  5V,  5T)  may  be 
combined  to  form  the  Lagrangian  and  collected  on  terms  with  common  variational 
displacement  factors. 


J^‘8Ldt  +  6L,|;|  =0  (3-30) 

^^lere  8L  =  8L5,  +8L5,  +8Lj,,  +5L,  and8L,  is  defined  at  die  time 

limits  t^  and  t, . 


Therefore, 
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=  .ffL  f£!^L+^'l+A„f5!t+^VA./^+^+2i^ 
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(3-33) 
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(3-35) 


The  boundaiy  term  of  the  Lagrangiaii  can  be  recast  in  the  following  general 
form  along  the  boundaiy. 


8Ls  =  J(N,8u,  +N^8u,  +Q,8w  +  M,8v,  +M„8\|;,)dS 


(3-36) 
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Figure  9:  Orientation  of  n  and  s  axes  along  the  plate  boundary. 

In  order  to  satisfy  Hamilton's  principle,  one  of  the  products  of  each  term  must 
be  specified  along  the  boundary  over  the  entire  time  interval.  This  leads  to  the 
following  conditions  \s4uch  must  be  specified  on  the  plate  boundaries.  The 
orientations  of  n  and  s  along  the  boundary  are  shown  in  Figure  9. 

N.oru,,  N^oru,,  Q,orw,  M,orv„  M^ory,  (3-37) 

The  final  term  of  the  Lagrangian  variation  (51^)  is  defined  only  at  the  limits  of 
the  time  integral.  In  order  for  this  term  to  satisfy  Hamilton's  principle,  the  generalized 
displacements  must  be  specified  at  the  time  interval  eni^ints.  This  condition  is 
satisfied  in  the  finite  element  method  by  discretizing  the  time  interval  and  treating 
each  time  step  in  a  semi-static  sense  as  explained  in  chiller  IV.  The  genmlized 
displacements  are  then  found  through  a  static  analysis  at  the  time  limits.  This  term  is 
shown  below  as: 


39 


dv  Sv, 

IT 


Sw 

8v  +  p,— 6w 

'  at 


dA 


(3-38) 


The  combination  of  all  these  variational  terms  yields  the  equation  form  given 
in  equation  3-39.  The  first  five  terms  in  the  area  integral  (first  line  below)  each 
include  tbe  integral  of  a  product  of  a  term  in  paiantheses  and  a  generalized  variational 
displacement  In  order  to  satisfy  Hamilton's  principle,  each  of  the  terms  in 
paiantheses  must  be  equal  to  zero  since  the  variational  displacements  are  aibitrary. 
This  (Hocess  generates  five  governing  equations  of  motion.  The  second  line  in 
equation  3-39  defines  the  boundaiy  and  initial  conditions. 


)8«  +  (  )5v  +  (  )8w  +  (  )8v,+(  )8\|;y]dA 
+J  [N.8u.  +  N.8u.  +  Q.8w  +  M,8v.  +  M,8\|;,  ]ds|  +  8L,  IJ;  =  0 


(3-39) 


By  iq)plying  the  defmition  of  strain  in  equation  3-2  and  the  load  resultant 
definitions  from  equations  2-8  and  2-9,  the  equations  of  motion  can  be  displayed,  in  a 
shorter  notation,  in  terms  of  load  resultants.  The  variational  displacement  before  the 
equations  below  describe  the  paratheses  location  of  the  corresponding  equation  in 
equation  3-39. 
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These  five  equations  are  the  governing  equations  for  all  plates  based  on  small 
deflection  theory  including  shear  deformation  and  rotary  inertia.  For  the  flnite 
element  method  presented  in  the  next  chiqrter,  a  displacement  formulation  is  required. 
The  equations  of  motion  can  be  cast  in  matrix  form  in  terms  of  displacements  for 
easier  conversion  to  the  finite  element  method.  This  leads  to  an  equation  of  the 
following  form: 


[M]{A}+[K]{A}  =  {f}  (3-41) 

wtere: 

{A}  =  {u,v,w,\|;,,M/yr 

{A}  =  {u.v,w,iir,,(irJ^ 


{f}  =  {0.0.q,0,0r 


(3-42) 


41 

The  mass  matrix  terms  are  shown  below: 

M„  =  M22  =  Mjj  =  p, 

M,4  =  Mjj  =  Ml,  =  Mj2  ^  P2  (3-43) 

M44  =  Mjj  =  p3 

All  other  M«  =  0 


The  stiffness  matrix  terms  are  written  as  differential  operators  on  the  vector 
{A}.  Note  that  this  matrix  is  symmetric  and  only  the  ui^)er  half  terms  are  indicated 
below. 


^11  =^1I^  +  3A,3^  +  A3j^ 

*^I2  =  A.,3^  +  (Au  +A33)^  +  A23^ 

K,3=0 

”  ®n^  +  2B,3^+B33^ 

^15  =  Bi3^(Bi2  +Bj3)^  +  B23^ 

K  —A  X--I-A  .^4.A 
^22  ~  ^33^^"23^y^"22  ^ 

K33=0 

^24  “  Bi3  (Bi2  B33  B23 

^25  =  B33  ^  +  2B23  ^  +  Bjj  ^ 

K33  =  A  JJ  ^  +  2  A4J  ^  +  A44 

^34  ~  •^35^'*'^43'^ 

K35  =  A4j4-  +  A44^ 

K44  =  D„^  +  2D,3^  +  D33^  +  A55 

^45  =  B,3'^  +  (D,2  +D33)-^+D23'^  + A4J 
Kss  =  D33  ^  +  2D23  ^  +  D22  ^  +  A44 


(3-44) 
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With  the  FSDT  presented  in  this  ch^iter,  the  ground  woit  is  set  to  implement 
this  theory  into  finite  elemmn  code.  The  next  chi4)ter  shows  how  the  above  matrix 
equation  is  discretized  into  finite  elements  and  tqpplied  to  the  computer  program. 


CHAPTER  rV 


FINITE  ELEMENT  FORMULATION 


This  chapter  describes  the  finite  element  method  (FEM)  formulation  used  in 
the  computer  program.  Most  of  the  FSDT  development  is  [Hesented  in  die  last 
chapter.  The  subsequent  FEM  formulation  follows  J.  N.  Redcfy  [25].  The  details 
shown  below  describe  how  the  fvevious  governing  equations  are  discretized  into 
elements.  It  is  assumed  that  the  reader  has  some  knowledge  of  the  finite  elemmt 
method  including  interpolation  functions. 

Generalized  Displacements  and  Interpolation  Functions 

For  the  finite  element  model,  the  domain  of  the  ]date  mid-plane  is  denoted  as 
R  and  is  divided  into  a  finite  number  of  elements  i^iose  domains  are  R,  (where  e  = 
1,2,3,...  number  of  elements).  For  each  element  domain,  R,  ,  the  generalized 
displacements  are  defined  by  use  of  an  interpolation  function  (where  i  =  1, 2, 3, ..., 
n).  The  interpolation  function  is  the  same  for  all  five  generalized  displacements 
^^ch  are  described 

»  =  v  =  Jvj^i  w  =  Jwi^i 

i-i  1-1  i-i 

•  a  ' 

1-1  1-1 
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n  is  the  number  of  nodes  per  element  and  the  interpolation  f  motion  depends  on  the 
type  of  element  used  in  the  analysis. 

Development  of  Element  Mam  and  Stiffness  Matrices 

The  equations  of  motion  derived  for  the  FSDT  can  be  aj^lied  to  each  element 
by  substituting  the  equations  4-1  and  into  equations  3-41  to  3-44.  The  simplified  form 
of  the  differential  equation  4-2  remains  the  same,  but  the  size  and  formulation  of  the 
matrices  and  vectors  changes.  The  mass  and  stiffness  matrices  are  square  symmetric 
matrices  of  the  order  five  times  the  number  of  nodes  per  element.  The  acceleration, 
displacement  and  force  vectors  are  of  the  same  order. 


[M']{4*}+[K’]{a'}  =  {f*} 


(4-2) 


where  the  matrices  and  vectors  are  defined  as: 


P,[S] 

0 

0 

P2[S] 
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P.[S] 
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0 
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e 
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P3[S]. 

(4-3) 


45 


[K-]  = 


K"- 

K** 

• 

f  t  ^ 

{F*} 

K“ 

K“ 

{F^} 

K“ 

K” 

K”‘ 

{f*}  =  . 

{FM 

K** 

K« 

{F^} 

K*' 

K« 

K” 

.{F’}. 

(-M) 


The  values  for  the  elements  of  the  stif&ess  matrix  are  given  below: 

[K"]  =  A„[S“]+A„([S->]+[S>'])+A„[S”] 
[k"]=a.4s*»]+a„[s”]+Aj,[s”]+A3,[s’-] 

[K->]  =  [0] 

[k«]=b,.[s“]+b„([s->]+[s>*])+b^[s»] 

[k‘’]=b,4s»]+b,3[s“]+b^[s»']+b„[s”] 

[K=']=[K“f  =  A,3[S’“]+A,3[S"]+A3,[S”’]+A„[S’>] 
[k=  ]  =  A33[s>’  ]+ A3,([s"'  ]+ [S’*  ])  +  A„[S”  ] 

[k’>]=[0] 

[k“  ]  =  B,3  [S’*  ]  +  B„[S“  ] + B3,[S”  ]  +  B„[S’’  ] 
[K“]=B33[S>’]  +  B33([S‘’]  +  [S>*])  +  B„[S-] 

[k"]=[o] 

[k”]=[o] 

[K”]  =  A,,[s>’]+A«([S’’]+[S’*])+A„[S“] 

[K“]  =  A„[S’»]+A„[S*»] 


(4-5) 
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[K“]  =  A„[S’«]  +  A„[S*>] 

[K"  ]  =  [K”  f  =  B„  [S“  ]  +  B.,([S-  ] + [S'  ]) + B„[S”'  ] 

[k'  ]  =  [k**  ]'  =  B„  [S'  ]  +  B„[S-  ]  +  B^IS”  ] + B„[S’'  ] 
[K'].[K"f  =A„[S']+A„[S"] 

[K']  =  D„[S“]  +  D,.([S']+[s>-])+D3,[S”]  +  A„[S] 

[K«]  =  D„[S'  ]  +  D„[S'  ] + D„[S'  ]+ D^IS'  ]  +  A„[S] 

[K”  ]  =  [K‘>f  =  B,4s-  ] + B„[S"  ] + B^[S”'  ] + B„[S'  ] 

[K“]  =  [K“]"=B^[S']+Ba([S']+[S'])+B„[S”] 

[k”]=[k“]'  =  a«[s“>]+a„[s'>*] 

[K“  ]  =  [K«  f  =  D„  [S’-  ]  +  D.,[S”  ] + D„  [S'  ] + D„[S'  ] + A„  [S] 
[K“].D„[S”]+D„([S']+[S>-])+D„[S”]+A«(S] 


The  values  of  the  sub-matrices  are  area  integrals  of  interpolation  function 
values  and  derivatives.  These  sub-matrices  are  of  order  n  x  n,  and  are  evaluated  by 
the  following  definitions  (i,  j  = 


S*i=J^  ♦i^jdxdy 
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Note  that  the  and  [k"  ]  sub-matrices  are  not  symmetric.  However,  since  diese 
matrices  have  the  (voperties  [S*^’]  =  [S’*]  and  [k“]  =  [k”]^,  the  total  element 
matrices  [m‘  ]  and  [k‘  ]  are  symmetric. 

The  sub-vectors  of  the  force  vector  are  of  order  n  and  are  ^ven  by  their 
components  (i  =  l,2,3,.,.4i): 

F“  =  f  f.4idxdy  +  Pj ,  a  =  1.2,3,4,5  (4-7) 

For  uniform  transverse  pressure  as  used  in  this  program,  fj  =  q  (the  transverse  {vessure 
value),  all  other  f,  0,  and  Pj  are  the  nodal  contributions  of  the  boundary  force 

conditions  along  the  plate  boundaries. 

Element  Types 


4— Node  8— Node  9-Node 

4  3473473 


1  2  15  2  15  2 

Figure  10:  Element  types  and  nodal  point  numbers. 


The  program  allows  the  user  to  choose  from  the  three  element  types  shown  in 
Figure  10:  four-noded  linear  quadrilateral,  eight-noded  quadratic  quadrilateral,  and 


nine-noded  quadratic  quadrilateral  elements.  The  four,  ei^t,  and  nine-noded 
elements  fxoduce  element  matrices  of 20x20, 40x40,  and  45x45,  respectively. 

The  interpolation  fimctions  used  for  these  elements  are  isoparametric  and 
belong  to  the  Lagrange  fiunily.  The  terms  of  the  [S^'i]  matrices  are  found  using  Gauss- 
Legendre  quadrature  numerical  integration.  Full-integration  is  used  for  all  stiffness 
terms  except  for  those  terms  involving  traverse  shear  coefficients  (A44,  A^j,  A,,)  in 
\^ch  a  reduced-integration  scheme  is  used.  The  reduced-integration  is  performed  to 
prevent  shear-locldng  effects.  J.  N.  Re(&]y  presents  a  more  complete  treatment  -f 
these  subjects  [25]. 

Finite  Element  Procedure 

For  static  problems,  the  differential  equation  becomes  a  much  more  simple 
series  of  linear  equations  because  the  acceleration  vector  is  zero.  In  this  case,  the 
following  equation  applies: 

[k']{a‘}  =  {f}  (4-8) 

In  this  instance,  the  element  matrices  are  assembled  globally  into  a  banded 
matrix.  Boundary  conditions  are  then  applied  to  the  global  equation  and  the  equation 
is  solved.  The  resulting  vector  gives  the  generalized  displacement  values  in  the 
following  order: 
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These  displacement  values  are  used  to  calculate  strains  and  stresses  at  tl^  Gaussian 
points  of  each  element  by  l^)plying  the  definitions  of  strain  and  constitutive 
relationships  equation  2-3  fiom  chqjter  n.  The  following  stresses  are  calculated  at 
each  lamina  interface  throug^i  the  [date  thickness  at  the  Gaussian  coordinates: 

Time-Dependent  Formulation 

Time-dependent  problems  require  some  extra  steps  in  the  solution  because 
they  require  the  solution  of  a  second-order  differential  equation.  This  method  follows 
that  presented  by  Reddy  [24^5].  The  {xocess  requires  transforming  the  equation 
shown  below  into  a  series  of  linear  equations. 

[M]{a}+[K]{a}.{f}  (4-10) 

One  way  of  tiansfonnii^  this  equation  into  a  solvable  form  is  by  using  discrete 
time  steps  in  the  analysis.  The  method  used  in  the  jxogram  is  known  as  the  Newmaric 
integration  scheme.  A  time  step  (At)  is  chosen  which  yields  a  stable  and  accurate 
solution  as  described  later.  For  each  time  step  equation  4-9  can  be  expressed  in  the 
following  general  discretized  form; 

iK]{A}„,={F}  (4-11) 

>^re  the  components  are  defined  as: 


[K]  =  [K]+a.[M] 
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{f}  =  {F}„,  +[M](a.{A}.  +.,{a}_  +ajA}  j 


(4-12) 


a,  =  aoAt  a,  - 1 

1  0  2  2p 


The  values  for  {A}^.{A}^.and{A}^arc  the  initial  velocities  and  wxelerations 

sui^lied  by  the  program  user.  Once  the  displacement  vector  at  the  new  time  step 
{A}^,  is  found  by  solving  equation  4-11,  the  new  accelerations  and  velocities  may 
calculated  using  the  following  equations: 


The  a  and  p  variables  are  parameters  used  in  the  Newmaric  scheme  to  create  a 
stable  and  accurate  integration  solution.  There  are  two  widely  used  pairs  for  these 
parameters  which  guarantee  stability  in  the  analysis. 

Linear  acceleration  method:  a  =  1/2,  p  =  1/6  (4-14) 

Constant  acceleration  method:  a  =  1/2,  P  =  1/4 

These  two  pairs  guarantee  stability  in  the  time  integration  scheme,  but  they  do 
not  necessarily  [Hovide  accuracy.  In  order  to  obtain  accurate  results,  the  time  step 
must  be  chosen  ai^)ropriately.  In  most  cases,  shorter  time  steps  yield  more  accurate 
approximations  than  longer  ones.  The  following  formula  provides  one  way  to  choose 
an  adequate  time  step. 
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(4-15) 

In  this  fonnula,  d  is  the  miniimiin  distance  between  any  two  nodal  points,  p,  is  the 
transverse  inertial  term  defined  in  equation  2-12,  and  D  is  the  lesser  of  the  two 
bending  stiffness  terms  D,|  and  Dj,  defined  in  equation  2-10.  Care  should  be  given  in 
choosing  too  small  a  time  step  because  of  the  computational  cost  of  that  choice. 

With  these  additional  steps,  the  procedure  presented  for  the  static  case  is 
followed.  This  procedure  requires  solving  a  set  of  linear  equations  for  the 
displacement  vector  and  computing  then  new  velocity  and  acceleration  vectors  for 
each  time  step.  For  this  reason,  the  computing  expense  can  become  quite  high. 

Boundary  Conditions 

Equation  3-37  dictates  the  boundary  condition  pairs  v^ch  must  be  specified 
along  the  plate  edges.  When  choosing  boundary  conditions  for  a  plate  problem,  one 
of  each  pair  must  be  specified  at  each  nodal  point 

N,  oru,,  N.oru,,  Q,  orw,  M,  orVa.  M^or\j/,  (4-16) 

The  following  list  describes  the  ^licable  force  and  displacement  values  for 
some  commonly  used  boundary  conditions. 

1.  Simply-Supported  Edge 

N.=0,  N^=0,  w  =  0,  M.=0.  M/.=0  (4-17) 

2.  Hinged  Edge  -  Free  in  the  normal  direction 
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N,  =0,  u,  =0,  w  =  0,  M,  =0.  M/,  =  0  (4-18) 

3.  Hinged  Edge  -  Free  in  the  tangential  direction 

u,  =  0,  =0,  w  =  0,  M,  =  0,  V,  =  0  (4-19) 

4.  Clamped  Edge 

u,  =  0,  u,  =  0,  w  =  0.  V,  =  0,  v,  =  0  (4-20) 

5.  Free  Edge 

N.=0.  N„=0.  Q.=0.  M.=0.  M„=0  (4-21) 

6.  Line  of  Symmetry  (for  symmetrical  finite  element  ixoblems) 

u.=0.  N^=0,  Q.=0,  v.=0.  M^=0  (4-22) 


Program  Information 

For  implementation  of  the  displacement  and  force  boundary  conditions  in 
COMPLATE,  the  following  information  is  necessary.  By  default,  all  generalized 
displacements  are  assumed  to  be  free  to  move  and  all  forces  are  assumed  to  be  zero. 
The  displacement  and  force  vectors  each  have  five  times  the  number  of  no(fes  entries. 
The  ordering  of  the  displacement  and  force  vectors  are  shown  below: 

Displacement  Vector -►{{u,v,w,\|i,,\|^y}“^,{u,v,w,\|;,,\|/y}“^,...} 

r  1  (4-23) 

Force  Vector 


To  change  a  boundary  condition  from  the  default  conditions,  the  user  must 
siq^ly  the  position  of  the  displacement  or  force  in  its  respective  vector  and  the  value 
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of  that  condition.  For  example,  to  specify  a  displacement  of  10  in  the  y-direction  on 
nocfe  3,  the  user  would  give  12  for  its  position  in  the  vector  and  10  for  its  value. 

Summary 

The  basic  theory  behind  the  development  of  COMPLATE  has  been  {xesented 
in  the  last  two  chapters.  This  {mrgram  is  written  in  the  FORTRAN-77  standard  for  use 
on  microcomputers  although  the  code  is  generic  and  may  be  used  on  larger  systems. 
Appendices  B  and  C  provide  more  information  on  the  program  including  user 
information,  sample  program  input  and  output  data  files,  and  documented  source 
code. 


CHAPTERV 


NUMERIC AL  EXAMPLES 


Overview 

In  this  chapter  several  numerical  examples  are  presented  to  validate  the 
computer  code.  Although  the  computer  program  can  analyze  more  comi^icated 
composite  plate  shapes  and  laminate  lay-ups,  the  following  cases  are  chosen  because 
other  solution  methods  are  available  for  these  types  of  composite  plates  and  loadings. 
Two  static  and  one  dynamic  plate  problems  are  considered.  The  static  problems  are 
compared  against  analytical  solutions  using  the  Navier  series.  The  dynamic  problem 
is  compared  against  a  solution  given  in  the  literature. 

The  Navier  series  solution  method  is  i^esented  in  ^^)endix  A  for  the  two 
plate  {xoblems  in  this  chapter.  The  Navier  series  is  used  to  find  an  exact  solution  to 
both  the  classical  plate  theory  (CPT)  and  first-order  shear  deformation  theory  (FSDT) 
plate  equations.  The  CPT  solution  demonstrates  the  importance  of  shear  deformation 
in  the  aruilysis. 

The  plate  used  in  all  three  cases  is  square  with  sides  of  length,  a,  as  shown  in 
figure  11.  The  loading  condition  for  all  three  cases  is  a  uniform  transverse  pressure 
aj^lied  over  the  area  of  the  plate  surface.  The  laminate  stacking  sequence  and  the 
boundary  conditions  change  for  the  three  cases  in  order  to  utilize  more  of  the 
computer  code. 
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Figure  11:  Plate  dimensions  and  loading  for  Cases  1, 2,  and  3. 


For  the  two  static  cases.  Cases  1  and  2,  the  following  material  properties  are 
used  v^ch  are  representative  of  composite  materials.  These  particular  properties 
were  introduced  in  1969  [IS]  and  have  been  used  by  many  other  author's  since. 


E,  =25x10 

G,2  =  Gi3  =  O.Sx  10‘psi 

v,j  =  0.25 


Ej  =  1 X 10*  psi 
Gjj  =  0.2x10*  psi 


(5-1) 


Case  1:  Symmetric  Specially-Ortfaotropic  Square  Plate  under  Uniform 
Transverse  Pressure  (Simply-Supported) 

For  this  case,  a  symmetric  specially-orthotropic  laminate  is  chosen  for  its 
laminate  property  features.  The  definition  of  such  a  laminate  is  one  >riiich  exhibits 
the  following  laminate  properties. 
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A, 3  —  A23  =  A4J  —  —  D|3  —  D23  —  0  (5-2) 

This  type  of  laminate  behaves  like  an  oithotroinc  material  for  gmnal  response 
purposes  because  the  laminate  stiffitess  matrix  contains  the  same  zero-valued  terms  as 
an  orthotroiHC  material.  All  the  bending-twisting  coupling  terms  vanisL  For  the 
calculation  of  interlaminar  stresses,  the  lamina  constitutive  relationships  are  required 
Therefore,  a  |xogram  ^ch  is  able  to  analyze  ordiotropic  plates  is  insufficient 
Because  the  B  matrix  is  zero  and  there  is  no  iq^ed  in-plane  loading,  in-plane 
displacements  decouple  from  the  transverse  deflection,  w,  and  the  rotations, 
M;,  and  Vy .  This  results  in  no  in-plane  displacements  of  the  laminate  mid-plane  (u  =  v 

=  0). 

For  this  case,  the  laminate  stacking  sequence  is  [0/90/0]  vdiere  all  three  layers 
have  equal  thicknesses,  luunely  h/3.  This  type  of  lay-up  is  called  a  cross-ply  because 
the  laminate  is  constructed  of  only  0  and  90  layers. 

This  plate  is  assumed  to  have  simply-supported  boundary  conditions  along  all 
the  edges.  This  results  in  the  following  boundary  conditions  (ignoring  in-plane 
boundary  conditions). 

w(0,y)  =  w(a,y)  =  w(x,0)  =  w(x,b)  =  0 

V,(x,0)  =  v,(x,b)  =  Vy  (0,y)  =  Vy(a,y)  =  0  (5-3) 

M,(0,y)  =  M,(a,y)  =  My(x,0)  =  My(x,b)  =  0 


This  case  is  used  to  evaluate  the  effectiveness  of  each  element  type  compared 
to  the  exact  analytical  solution  so  all  three  element  types  are  used  in  the  analysis. 
Also  the  effect  of  shear  deformation  is  observed  by  using  three  length  to  thickness 
ratios  (a/h  =  100,10,4). 
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Because  of  this  problem's  symmetry,  only  one  quarter  of  the  plate  is  needed  for 
the  analysis.  This  results  in  the  meshes  shown  in  figure  12. 


8  and  9-Noded  4-Noded 

Element  Meshes  Elonent  Mesh 

(A.B) 


a 


Figure  12:  Geometry  of  FEM  mesh  for  Cases  1  and  3. 
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The  results  of  this  analysis  are  shown  in  Table  1.  The  deflection  and  stresses 
are  non-dimensionalized  by  the  definitions  given  in  equation  5-4.  The  stresses  shown 
in  the  table  are  calculated  at  Gaussian  points  with  coordinates  defined  by  A  and  B  in 
equation  5-5. 
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■'b  =  ''b(B,A) - in  layer  3  (0®) 

qo® 


=  Ty^(A,B)-—  in  layer  2  (90®) 
Qo® 


A=  0.55283  X  a  (for  2x2QX  0.56250  x  a  (for  4x4L) 

B=  0.94717  X  a  (for  2x2Q),  0.93750  x  a  (for  4x4L)  (5-5) 


The  terms  in  the  "Method"  colunm  of  Table  1  refer  to  following  ^  methods  of 
analysis: 


COMPLATE  Results: 

2x2Q9-  4  element  mesh  using  9-noded  quadratic  elements. 

2x2Q8-  4  element  mesh  using  8-noded  quadratic  elements. 

4x4L4- 16  element  mesh  using  4-noded  linear  elements. 

Other  method  results: 

CPT-  Navier  series  solution  using  classical  plate  theory.  Appendix  A. 

FSDT-  Navier  series  solution  using  first-order  shear  deformation  theory.  Appendix  A 
HSDT-  Center  deflection  solution  by  Reddy  using  higher-order  shear  (^formation 
theory[22]. 

Navier  series  solutions  utilized  49  terms  (m,n  =  1...49)  -  see  Appendix  A. 


The  results  in  Table  1  show  that  COMPLATE  produces  results  close  to  those 
of  the  Navier  solution  using  FSDT  for  coarse  meshes  using  both  8  and  9-noded 
quadratic  elements.  The  4-noded  linear  element  appears  to  be  less  accurate.  Also 
note  that  the  non-dimensionalized  deflection,  w,  is  much  higher  for  the  cases  a/h=10 
and  4  than  predicted  by  CPT. 
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Table  1;  Comparison  of  maximum  deflection  and  stresses  for  Case  1  [0/90/0], 


— 

Method 

W 

^yz 

CPT 

06660 

0.7908 

0.1935 

-0  03944 

0 

0 

2x2Q9 

0.6704 

0.7931 

0.1949 

-0.03861 

-0.7036 

-0.2049 

2x2Q8 

0.6707 

0.7933 

0.2015 

-0.03757 

-0.7118 

-0.1999 

a/h=100 

4x4L4 

0.6660 

0.1152 

0.1960 

-0.03631 

-0.6968 

-0.1967 

HSDT[22] 

0.6705 

FSDT 

0  6697 

0.7905 

0.1948 

-0,03929 

-0.7020 

-0.2020 

2x2Q9 

1.0234 

0.7577 

0,3076 

-0.04555 

-0.6865 

-0.2336 

2x2Q8 

1.0211 

0.7575 

0.3078 

-0.04516 

-0.6872 

-0.2340 

a/h=10 

4x4L4 

1.0276 

0.7386 

0.3096 

-0.04332 

-0.6758 

-0.2247 

HSDT[22] 

1.0900 

FSDT 

1.0219 

0.7556 

0.3066 

-0.04657 

-0.6823 

-0.2294 

2x2Q9 

2.6630 

0.6419 

0.6513 

-0.06588 

-0.6143 

-0,3275 

2x2Q8 

2.6559 

0.6419 

0.6513 

-0.06544 

-0.6143 

-0.3278 

a/h=4 

4x4L4 

2.7025 

0.6238 

0.6494 

-0.06300 

-0.6020 

-0.3198 

HSDT[22] 

2.9091 

FSDT 

2.6595 

0.6408 

0.6494 

-0.06725 

-0.6104 

-0.3230 
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Case  2:  Angle-Ply  Square  Plate  under  Uniform  Transverse  Pressure  (Hinged) 

For  this  case,  an  angle-ply  laminate  is  chosen  for  its  laminate  property  features. 
The  definition  of  such  a  laminate  is  one  with  the  stacking  sequence  [+6/-0].  where  n  is 
some  integral  multiple.  This  type  of  laminate  has  the  following  laminate  property 
simplifications: 

A|3  =  A23  =  A^j  =  B,j  =  B|2  =  B22  =  B33  =  Dj3  =  D23  =  0  (5-6) 

The  only  bending-twisting  coupling  terms  appear  in  the  B  matrix.  This  case  is 
able  to  test  the  program's  coiq)ling  effect  calculations. 

For  this  case,  the  laminate  stacking  sequence  is  [4S/-4S]2  where  all  four  layers 
have  equal  thicknesses,  namely  h/4.  This  plate  is  assumed  to  have  hinged  edges  with 
freedom  in  the  tangential  direction  along  the  plate  boundaries.  This  results  in  the 
following  boundaiy  conditions: 


w(0,y)  =  w(a,y)  =  w(x,0)  =  w(x,a)  =  0 
u(0,y)  =  u(a,y)  =  v(x,0)  =  v(x,a)  =  0 

V,(x,0)  =  M/,(x,a)  =  M/y(0,y)  =  v,  (a,y)  =  0  (5-7) 

My  (x,0)  =  My  (x,a)  =  M,(0,y)  =  M,(a,y)  =  0 
N.(x,0)  =  N.(x,a)  =  Ny(0,y)=  Ny  (a,y)  =  0 


Since  all  three  element  types  were  evaluated  in  the  preceding  case,  only  the  9- 
noded  quadratic  element  is  used  for  this  case.  B  cause  the  laminate  stacking 
sequence  is  not  symmetric  about  its  mid-plane,  the  same  plate  symmetry  as  in  case  1 
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does  not  exist  for  this  case.  Therefore  a  full  4x4  mesh  of  the  entire  plate  is  used  as 
shown  in  Figure  13. 
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Figure  13:  Geometry  of  FEM  mesh  for  Case  2. 

The  results  of  this  analysis  are  shown  in  Table  2.  The  deflection  and  stresses 
are  non-dimensionalized  by  the  definitions  given  in  equation  5-8.  The  stresses  shown 
in  the  table  are  calculated  at  the  Gaussian  points  coordinates  (tefined  by  A. 


(5-8) 
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A=  0.55283  X  a  (5-9) 

4x4(^  refers  to  results  from  CX)MPLATE.  The  other  methods  are  defined  in 
the  same  way  as  in  case  1.  In  addition  to  tabular  results,  the  stresses  are  plotted 
through  the  thickness  to  show  a  representation  of  interlaminar  stresses  in  Figures  14  to 
16. 

The  results  in  Table  2  show  that  the  computer  im>gram  provides  adequate 
accuracy  compared  to  the  Navier  solution  of  FSDT.  Also  note  that  the  non- 
dimensionalized  deflection,  w,  is  much  higher  for  the  cases  a/h=10  and  4  than 
predicted  by  CPT. 
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Table  2:  Comparison  of  maximum  deflection  and  stresses  for  Case  2  [±45]2 . 


Method 

w 

o. 

CPT 

0.4408 

0.2040 

0.2040 

-0.1774 

0 

0 

a/h“100 

4x4Q9 

0.4439 

0.2041 

0.2041 

-0.1776 

-0.01783 

-0.01783 

FSDT 

0.4433 

0.2039 

0.2039 

-0.1773 

-0.01789 

-0  01789 

«/h=10 

4x4Q9 

0.6925 

0.2011 

0.2011 

-0.1751 

-0.01773 

-0.01773 

FSDT 

0.6917 

0.2007 

0.2007 

-0.1746 

-0.01788 

-0.01788 

a/h=4 

4x4Q9 

2.0186 

0.1986 

0.1986 

-0.1731 

-0.01763 

-0.01763 

FSDT 

2.0164 

0.1977 

0.1977 

-0.1722 

-0.01784 

-0.01784 

«»x 


Figure  14:  Normal  stress,  a.,  through  the  plate  thickness 
near  the  plate  center  for  Case  2  [±45]2 . 
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Case  3:  Orthotropic  Square  Plate  under  Suddenly  Applied  Uniform  Transverse 
Pressure  (Simply-Supported) 

In  order  to  validate  the  dynamic  analysis  routines,  a  simple  orthotroiHC 
laminate  is  chosen  for  this  case.  Results  presented  by  Reddy  [23]  are  used  in 
comparisoa  The  material  {voperties  and  time  parameters  used  in  this  analysis  are 
also  taken  from  [23].  The  following  material  properties  are  used 


and  the  time  panuneters  for  the  Newmark  scheme  (constant-average  acceleration)  are: 

At  =  5  psec,  a  =  0.5,  P=0.25  (5-11) 

Since  only  one  layer  exists,  the  stacking  sequence  is  [0].  The  plate  edges  are 
simply-suj^rted  and  a  uniform  pressure  is  suddenly  iq>plied  at  the  initial  time  step 
wdiich  remains  steacfy  throughout  the  analysis.  The  boundary  conditions  are  the  same 
as  in  Case  1.  Symmetry  is  again  utilized  and  the  mesh  is  the  same  as  the  one  used 
with  the  2x2Q9  elements  in  Case  1.  Since  a  direct  comparison  is  made  with  the 
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literature,  the  results  are  not  non^limensionalized  The  plate  is  square  with  tiie 
following  geometric  parameters; 

a  =  25  cna,  h  =  5  cm  (5-12) 

Table  3  shows  the  results  of  the  center  deflection  and  the  stress  at  the  Gaussian 
point  nearest  to  the  plate  center  (A,A)  over  time  compared  to  the  results  i^-esented  by 
Reddy.  Since  Reddy  used  a  similar  {vocedure,  the  results  coincide  exactly  except  for 
a  couple  of  what  a]q)ear  to  be  typognq^cal  errors  in  his  article.  This  correlation 
validates  the  (^mamic  routines  of  the  fvogram. 
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Table  3;  Dynamic  response,  deflection  (w)  and  nonnal  stress  (o,),  of  a  square 
orthotroinc  plate  under  pulse  uniform  transverse  pressure.  Case  3. 


time 

a«ec) 

w-10*  (cm) 

ReddYr231  2x2Q9 

Red(lvf231 

(N/Clll2) 

2x209 

10 

0.0079 

0.007963 

2.986 

2.986 

20 

0.0398 

0.03985 

24.64 

24  64 

40 

0.1939 

0.1939 

132.2 

132.2 

60 

0.4303 

0.4303 

282.1 

282.1 

80 

0.5531 

0.5531 

359.3 

359.3 

100 

0.5264 

0.5264 

349.7 

3497 

120 

0.3705 

0.3705 

345.4 

245.4 

140 

0.1779 

0.1779 

115.1 

115.1 

160 

0.0353 

0.03533 

22.0 

22.0 

180 

-0.0395 

-0.03946 

-20.97 

-20.97 

200 

0.1105 

0.1105 

73.61 

73.61 

220 

0.3296 

0.3296 

214.1 

214.1 

240 

0.4781 

0.4781 

316.8 

316.8 

260 

0.5548 

0.5548 

368.9 

368.9 

280 

0.4797 

0.4797 

314.5 

314.5 

300 

0,2006 

0.3006 

194.9 

194.9 

320 

0.0840 

0.08402 

59.38 

59.38 

340 

-0.0302 

-0.03020 

-18.53 

-18.53 

360 

0.0459 

0.04587 

28.57 

28.57 

CHAPTER  VI 


CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 

The  primaiy  objective  of  this  work  was  to  develop  a  program  to  analyze 
general  laminated  composite  plates  under  static  and  dynamic  conditions  using  first- 
order  shear  deformation  theory.  Several  composite  plate  cases  were  used  to  validate 
the  computer  code  by  comparing  the  results  against  other  solution  methods  and  results 
found  in  literature.  The  following  conclusions  may  be  drawn  from  the  results: 

1.  Based  on  the  three  cases,  the  inogram  aiq3ears  to  correctly  wppXy  the  first- 
order  shear  deformation  theory  of  composite  plates  to  the  finite  element  method. 

2.  Of  the  three  element  types  used  in  the  program,  the  8  and  9-noded  quadratic 
elements  give  the  best  results  with  neither  one  showing  clear  superiority.  The  4-noded 
linear  element  ai^xars  to  give  less  than  adequate  results  compared  to  the  two  previous 
element  types. 

3.  The  deflections  predicted  by  the  computer  program  show  that  shear 
deformation  effects  can  be  significant  beyond  the  usual  range  for  isotropic  materials 
(a/h  >  10).  This  is  a  clear  demonstration  for  the  need  of  a  shear  deformation  theory 
for  the  analysis  of  all  but  very  thin  composite  plates. 

4.  Several  authors  have  presented  a  complete  treatment  of  classical  plate 

theory  for  composite  plates.  Because  of  length  restrictions  in  published  articles, 

FSDT  is  often  {^resented  in  abbreviated  form.  Current  composite  textbooks  give  only 

68 


69 


a  cursoiy  treatment  of  the  first-order  shear  deformatioD  theory  and  i^glect  to  detail  the 
full  importance  of  shear  deformation  in  composite  plates.  This  worii  rqxesoits  a  full 
presentation  of  this  theory  combined  into  a  single  source. 

Recommendations 

Based  on  research  in  composite  plate  theories  and  the  development  of  the 
computer  {vogram,  the  following  modifications  and  future  research  are  anticipated: 

1.  This  (xogtam  should  be  used  in  advanced  composite  classes  to  give 
students  insight  on  the  importance  of  shear  deformation  in  composite  plates  and  to 
gain  a  better  understanding  of  the  finite  element  formulation  of  composite  plate 
theories. 

2.  This  program  could  be  used  in  conjunction  with  composite  testing  if  funds 
are  available  for  future  work.  The  fxogram  could  be  used  in  developing  the  size  and 
stacking  sequence  of  the  composite  test  specimen,  and  also  the  fvogram  results  could 
be  directly  compared  with  test  results. 

3.  Because  the  environmental  responses  of  composites  are  significant,  thermal 
and  hygral  (moisture)  effects  should  be  incorporated  in  future  revisions  of  this 
I»ogram. 

4.  The  program  la’ovides  easy  ^>plication  of  a  uniform  transverse  pressure. 
Other  types  of  transverse  loading  require  a  significant  amount  of  work  on  the  user's 
pert  A  subroutine  should  be  developed  to  ai^ly  more  general  types  of  loading 
functions,  q  =  q(x,y). 

5.  This  program  provides  a  good  starting  point  to  develop  a  program  based  on 
higher-order  shear  deformation  theories  of  composite  plates.  Because  of  the 
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additional  degrees  of  freedom  required  in  these  theories,  a  microcomputer-based 
ixogram  ^^)ears  to  be  prohibitive  at  fxesent  However,  a  program  can  be  developed 
on  a  larger  system  or  in  the  future  on  smaller  computers  as  memory  management 
problems  are  solved. 


APPENDIX  A 


NAVIER  SERIES  SOLUTIONS 


Navier  fvoposed  using  a  double  Fourier  series  to  solve  certain  differential 
equations.  This  type  of  series  was  applied  to  the  solution  of  plate  governing  equations 
with  particular  sh{q)es  and  boundary  conditions.  Because  this  solution  satisfies  the 
plate  governing  equations  and  meets  the  applicable  boundary  conditions,  it  is  viewed 
as  an  exact  solution  of  the  plate  theory  used  in  the  analysis.  Since  the  series  is 
infinite,  it  is  impossible  to  utilize  all  the  terms.  However  taking  enough  terms  in  the 
solution  (m,n  <  50),  the  solution  converges  to  the  significant  figures  presented  as 
results  in  chapter  V. 

For  this  analytical  solution,  general  transverse  pressure,  q(x,y),  on  a 
rectangular  plate  can  be  tqrplied  by  use  of  the  following  Navier  series: 


q(x,y)  =  2ZQ-s“» 


imcx  .  mcy 

- sm — - 

a  b 


(A-1) 


The  loading  terms  (Q,„  )  can  be  found  by  integrating  equation  A-1,  and  then 
solving  for  Q„ .  This  yields  the  following  result: 


Q™  =  TJn  Jn  - sm— ^dydx 

ab  •'0  •'0  a  b 


(A-2) 
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Since  the  cases  presented  in  the  results  utilize  a  uniformly  ap>lied  transverse 
pressure  (q(x,y)  =  qj,  it  is  convenient  to  defuie  the  loading  terms  for  this  condition. 
For  a  uniform  transverse  pessure,  these  terms  result  in  the  following  expression; 


^  q 

— -^(l  -  cosmx)(l  -  cosnji) 
mmt 


(A-3) 


CASE  1:  SYMMETRIC  SPECIALLY-ORTHOTROPIC  LAMINATES 


Symmetric  specially-orthotropic  plates  are  often  used  to  compare  results 
because  they  simplify  the  governing  plate  equations.  For  plates  with  simply-supwrted 
edges,  the  analysis  pocess  can  be  accomplished  by  use  of  a  Navier  solution.  This 
type  of  plate  is  defined  as  being  symmetric  with  respect  to  the  laminate  mid-plane  and 
whose  bending-twisting  coupling  terms  vanish.  This  translates  to  plate  lay-up  with 
the  following  material  poprty  simplifications: 

Aj3  =  A23  =  Bjj  =  Djj  =  D23  =  A4J  =  0  (A-4) 


Classical  Plate  Theory  Solution: 

The  CPT  solution  ignores  shear  deformation  and  is  therefore  indepndent  of 
thickness  effects.  This  solution  is  used  as  a  baseline  case  for  very  thin  plates. 
Whimey  [33]  provides  a  more  complete  treatment  of  this  method.  The  following 
boundary  conditions  are  used  for  the  CPT  solution  of  the  simply-supported  rectangular 
plate: 
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w(0,y)  =  w(a.y)  =  w(x,0)  =  w(x,b)  =  0 

^(x.O)  =  ^(x,b)  =  ^(O.y)  =  ^(a.y)  =  0 

OK  OK  dy  ay 

M.(O.y)  =  M,(a.y)  =  M,(x.O)  =  M,(x,b)  =  0 


(A-5) 


With  the  given  material  property  simplifications,  the  governing  CPT  equations 
reduce  to  the  following; 


n^4  V  u  ^’dK^dy^  ^  dK*  ^ 


(A-6) 


The  deflection  function  is  chosen  to  satisfy  the  plate  boundary  conditions  and 
the  simplified  CPT  governing  equations  and  is  given  as;; 


A  X  .  nrty 

w(x,y)  =  sm - sm-^ 

a  b 


(A-7) 


By  applying  this  function  to  the  CPT  equation,  the  following  relation  for  the 
constant  is  obtained  (with  R  =  a/b): 


W_  = 


a^  Q. 


7l" 

where  =  D,,m^ +2(D,2 +2Dj^)(mnR)^ +D22(nR)^ 


(A-8) 


The  in*plane  stresses  in  the  layer  are  defined  by  the  following  equations; 


of  (x.y..) = +q»>„=r= 

a  b 
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o‘*Hx.y,z)  =  +Q^’n"R^)sin-^sm^ 


(A-9) 


x<«(x.y.z)=-2li<5S^XZ 

X  Ml  Ml 


mnQ^  nwtx  mcy 

-^^-5-^cos— cosY 

••I  •-!  l-'m»  ® 


Note  that  CPT  assumes  that  the  transverse  shear  strains  are  negligible,  and 
therefore  the  transverse  shear  stresses  are  zero. 


First-Order  Shear  deformation  theory  Navier  solution: 


Assume  the  following  Navier  Series  displacements  functions; 

w(x,y)  =  sin— sin-^ 

a  b 


,  \  ^  v'  iT«  mrtx  .  miy 

V,(x,y)  =  X2^^««cos - sm-^ 

■-I  a-l  a  D 


(A-IO) 


Vy  (x,y)  =  2^2-^’’- 

m-I  Ml  a  D 


These  assumed  displacement  functions  meet  the  followii^  simply-suf^rted 
boundary  conditions  on  the  plate; 


w(0,y)  =  w(a,y)  =  w(x,0)  =  w(x,b)  =  0 

M/,(x,0)  =  \ir,(x,b)  =  \|/,(0,y)=H/,(a,y)  =  0  (A-11) 

M,(0,y)  =  M,(a,y)  =  My(x,0)  =  M,(x.b)  =  0 


These  assumed  displacemer  t  fields  can  be  substituted  into  the  governing 
FSDT  equations  (3-40  to  3-43).  Since  the  B  matrix  terms  are  zero,  the  in-plane 


displacements  (u,v)  become  uncoupled  from  the  other  generalized  displacements  (w, 
iM'y  )•  Under  static  conditions  with  only  transverse  fvessure  fq^lied,  equations  3-40 
to  3-43  reduce  to  the  following  equations; 


^11  Ljj  L,j  ® 

L,j  Ljj  Ljj  -  %„  •  = .  0  >  (A-12) 

.^13  ^23  Ljj^  .Qb». 


^11  “  U,,a,  -J-Djjpjl  +  Ajj 
^12  ~  (Ui2  ■*■^33)01^^^ 

^13  = 

L22  +  D22Pa  +  A44 

L33  ~  Ass®a  ■*‘'^*«Pb 


.  mji  .  _  nrc 

>raere  a„  =  —  and  B.  =  — 
a  b 


(A.13) 


The  solution  of  this  equation  yields  the  constants  (W„, vydiich 

can  be  substituted  into  equation  A-10  to  determine  the  deflection  and  rotations  of  any 
point  on  the  plate  mid-plane.  The  plate  stresses  are  defined  as; 

o^‘’(x,y,z)=  -z2£(Q[‘>'P.„a.  -»-Q{‘>%_P,)sina,xsinP.y 

■-1  a-l 

oi^Hx.y.z)  =  -z^£(Q<‘>'P,^a,  +Q^^'Fy^P.)sina.xsinp,y 
■-1  0-1 


^SHx,y,z)=  z^  +'Py„a,)cosa„xcosP,y 


(A-14) 
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«0  «0 


xS’(x.y.2)  =  +W^PJsma.xcosp,y 

■-1  a-l 

+QS>(4'._  +  W_a.)cosa.xsiiip.y 


TS’(!‘.y.z)=  i;i;<5S’(%«  +  W_p.)smo.xa)sP.y 

B-l  a-t 

+Q<‘’('Pa„  +W^a,)cosa,xsinp.y 


CASE  2:  ANGLE-PLY  LAMINATES  (±6], 

Angle-ply  laminated  plates  are  also  used  to  compare  results  because  they  allow 
the  governing  plate  equations  to  be  simplified,  but  still  include  bending-stretching 
coiq)ling.  With  certain  boundary  conditions,  the  analysis  process  does  not  require  a 
numerical  solution.  Angle-ply  laminates  are  defined  by  a  laminate  stacking  sequence 
with  fiber  orientations  alternating  between  +0  to  -6  among  consecutive  layers.  The 
code  for  such  laminates  is  [±6],  or  [+6/-6],  ,  where  n  is  some  integer  multiple.  This 
translates  to  plate  lay-iq)s  with  the  following  material  property  simplifications: 

A|3  =  A23  =  B,,  =  Bj2  =  B22  =  B33  =  D|3  =  D23  =  A^j  =  0  (A-15) 


For  this  analytical  solution,  general  transverse  pressure  on  a  square  plate  is 
applied  in  the  same  manner  as  given  in  equations  A-l  to  A-3. 

JliJL  m«v  nm/ 

q(x.y)=2;2;<3- 

■-I  ■-! 


.  mm.  .  nm 

sm - sm — - 

a  b 


(A-16) 
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Classical  Plate  Theory  Solntion: 

The  CPT  solution  ignores  shear  deformation  and  is  independent  of  the  plate 
thickness  effects.  This  solution  is  used  as  a  baseline  case  for  very  thin  plates. 
Whiti^  [33]  provides  a  more  complete  treatment  of  this  solution.  The  plate 
consi(tered  here  is  rectangular  with  hinged-edges  (fm  in  the  tangential  direction).  In 
this  case,  the  following  boundary  conditions  apply  along  the  edges: 


w(0,y)  =  w(a,y)  =  w(x,0)  =  w(x,b)  =  0 
u(0,y)  =  u(a,y)  =  v(x,0)  =  v(x,b)  =  0 


^(0.y)  =  ^(..y)  =  ^(x.0)  =  ^(x.b)  =  0 
ay  ay  OK  ax. 

(O.y )  =  (a.y )  =•  (x.O)  =  (x.b)  =  0 

M,(0,y)  =  M,(a,y)  =  My(x,0)  =  M,(x,b)  =  0 


(A-17) 


With  the  material  simplifications  of  angle-ply  laminates,  the  governing  CPT 
differential  equations  become: 


K 

K 

K 


11 

12 

13 


K.2 

_ 1 

u 

0 

K22 

K23 

V 

0 

K23 

K33. 

w 

.q. 

(A-18) 


where: 


V  _  A  jf-4.  A 
'^ll  -  ^11  ^^33  jyJ 

^12  ~  (Ai2  Ajj)-^ 


*^13  -  +023  ^ 


K  —  A  -^  +  A 

^22  “  ^33  a|i  +^22  ayj 


(A-19) 
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By  substituting  these  displacement  functions  into  the  CPT  equations,  the 
following  matrix  equation  results: 

L,.  L,2  L,3lfU„^  0  ■ 

Li2  L22  L23  ]v„  [  =  1  0  .  (A-21) 

.^13  ^23  L33_  Qb*. 

where: 

L„  =Ai,ai+A33p^ 

^12  =  (^12  ■^•^33)®bPb 
Li3  =  ”3B,3a*  P,  -  B23P, 

^22  ~  ■^^22Pb 

L23  “  “®13®ii.  “  3®23®«Pn 
L33  =  D„al  +2(D,2  +2D33)aiPi  +D22P: 


(A-22) 
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.  mn  .  -  iwc 
and  a,  =  —  and  p,  =  — . 

a  b 


The  solution  of  these  equations  yields  the  constants  (U^,V^,W^)  which 
can  be  substituted  into  equation  A-20  to  determine  the  deflection  and  rotations  of  any 
point  on  the  plate  mid-plane.  These  constants  are  also  utilized  to  calculate  the  plate 
stresses.  The  in-plane  stresses  in  the  layer  are  defined  by  the  following  equations; 

o«>(x,y,z)=  +5g>V_p.  -2zf5„W_o.p.)  cosa.xcosP.y 

mal  a>l 

+(^w„  (Q,'f’ai  +  Qil’Pi  )-  Q!J’(U_P.  +  V_a.  ))suia.xsinp.y] 


<Ti*>(x,y.x)  =  i;  t  [(QS>U_o.  +  Q»>V„P,  -  2zQ„W„a.p. )  cosa,xcosp,y 

oi«l  n«l 

+(xW_((3g'a’.+Qg>p;)-(5g>(U_p.+V„a.))sina.xsinp.y] 


xS’(x.y.x)  =  i;i;[(Qg'U_a.+Q0>V_p.-2zQ„W_a.p.)  cosa„xcosp,y 
01-1  1-1 

+(xW_(Q2>a=.  +Q2’PI)-Q2’(U_P.  +  v_a.))siiia.xsmp.y] 


(A-23) 


Again  according  to  CPT,  the  transverse  shear  stresses  are  zero. 

First-Order  Shear  Deformation  Theory  Solution 

The  plate  considered  here  is  rectangular  with  hinged-edges  (free  in  the 
tangential  direction).  In  this  case,  the  following  boundary  conditions  apply  along  the 
edges; 
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w(0,y)  =  w(a.y)  =  w(x,0)  =  w(x,b)  =  0 
u(0,y)  =  u(a.y)  =  v(x,0)  =  v(x,b)  =  0 

\l/,(x,0)  =  M/,(x,b)  =  M/,(0,y)  =  Vy(a.y)  =  0  (A-24) 

M,(0,y)  =  M,(a,y)  =  M,  (x.O)  =  My(x,b)  =  0 
(O.y)  =  (a,y )  =  (x.O)  =  (x.b)  =  0 


This  solution  requires  live  displacement  functions.  The  following  Navier 
series  functions  satisfy  the  FSDT  governing  plate  equations  in  chapter  III  and  the  plate 
boundary  conditions: 

u(x,y)  =  K 

a-t  a-t  ®  D 

.  .  mxx  .  njty 

v(x,y)  =  2.2-^-.  - sm-j^ 

&  D 


w(x,y)  =  sin^^sin-^ 

a  b 


Wtml  a«l 


(A-25) 


/  ^  vv'iTi  “’’K  •  nxy 

V.(x,y)  =  2.2.^*-  - sm-^ 

a  b 


Vy  (x,y)  =  SZ^a.  sin-5^cos^ 

BMl  0.1  a  D 


Substituting  these  displacement  fields  into  the  governing  FSDT  equations  340 


to  343  yields: 
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Lm 

L.2 

L.3 

L,4 

L.,' 

0 

L.2 

L22 

L23 

^24 

^25 

0 

^13 

L23 

L33 

L34 

L33 

Q- 

L.4 

L24 

L34 

L44 

L43 

'5'*- 

0 

L,3 

L23 

L33 

^45 

1-53 . 

0 

(A-26) 


where  the  matrix  terms  are  defined  as: 
L|i  “ 

Li2  =  (Ai2  +  A33)<x^P, 

L,3  =  0 

Li4  =  2Bi3a^P, 

L22  =  Ajatti  +  Aap". 

L23=0 

L33  “  Ajjtt^  +  A44P2 
L34  “  Ajjtt^ 

L35  “  A44P, 

^44  ~  ^11®  m  ■•'^33P«  +Ajj 
=  (^12  ■*'I^33)®«iPii 
LjS  =  +D22PB  +  A^^ 


(A-27) 


and  a. 


The  solution  of  this  equation  yields  the  constants  (U„ ,  ,W„  ) 

^ch  can  be  substituted  into  equation  A'25  to  determine  the  deflection  and  rotations 
of  any  point  on  the  plate  mid-plane.  The  stresses  may  also  be  found  by  utilizing  these 
terms  are  defined  as  follows: 
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<j“>(x,y,i)  =  i;i;[(Q,‘f‘U_o.  +  Q!;'V_p.  +  zQSM^'.^P.  +  Vy.o.rtcosa.x 

m>l  a-l 

«»P.y-(5!;’(U_P.  +  V_a.)+2((5|;>T._a.  +  Q|,‘>'i'._p.)) 
sina.xsinP.y] 


o<«(x,y,z)  =  j;|;[{Q,TU_o.  +Qg'v_p.  +  2Q®(T._P.  +  4',_o.))«>sa.x 

m>l  >-I 

«>sp.y-fe‘’(U_P.  +  V_a.)  +  z(5,T'J'-_<».  +0S5”I'»_P.)) 

sina.xsin^.y] 


t!;'(x,y,2)  =  I;  j;[{Q<!'U_a.  +Q2*V_P.  +2Q<3‘'(4'.„P.  +  4',_a.))coso.x 
01-1  0-1 

a>sp.y-(Q2’(U_P.  +V_a.)+z{Qg>4'._o.  +  Qg>'F,_P.)) 
sina.xsinp.y] 


xS’(x,y,z)  =  J  +W„p.)sina.xcosP,y +  Q<$K'P»«o  +W„aJ 

0»-l  0-1 

cosa.xsinp.y] 


xi^’(x.y,z)=  +W^P,)sina,xcosp,y +  +W^a,) 

m-l  0-1 

cosa.xsinp.y] 


(A-28) 


APPENDIX  B 


COMPUTER  PROGRAM  INFORMATION 
COMPLATE  Program  Overview  and  Inatmctions 


Figure  17:  Defmition  of  coordinates  on  a  rectangular  plate. 

This  program  was  developed  to  analyze  general  laminated  composite  plates  of 
uniform  thickness,  h.  The  program  uses  first-order  shear  deformation  of  laminated 
orthotropic  plates  in  the  analysis.  The  results  from  the  program  include  the  laminate 
material  property  matrices,  [A],  [B],  [D],  the  resulting  mid-plane  displacements  and 
rotations  (u,  v,  w,  \\f^  \|iy)  at  plate  nodal  points,  the  laminate  strains  at  the  mid-plane 
and  on  the  laminate  sur&ces,  and  the  interlaminar  stresses  at  lamina  interfaces.  The 
following  paragraf^  describe  the  user  supplied  requirements  for  the  program.  Figure 
17  shows  the  coordinate  system  for  an  example  rectangular  plate. 
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1.  Program  Control  Parameters 

The  program  uses  several  parameters  to  control  the  scope  and  output  of  the 
program.  TITLE  is  a  one-line  output  heading  used  to  describe  the  laroblem  being 
analyzed.  The  next  two  parameters  dictate  the  element  type  used  in  the  jvogram.  lEL 
describes  the  order  of  the  element  (1  -  linear,  2  -  quadratic).  NPE  is  the  number  of 
nodes  per  element  (4  if  lEL  =  1,  8  or  9  if  lEL  =  2).  IMESH  is  the  parameter  used  to 
control  whether  a  rectangular  mesh  is  generated  (IMESH  =  1)  or  all  mesh  information 
is  input  for  a  general  sluq)e  mesh  (IMESH  =  0).  NPRNT  is  used  to  control  tl^  printing 
of  element  stiffness  matrices  and  force  vectors  (0  -  not  fainted,  1  -  printed).  ITEM 
indicates  whether  the  analysis  is  static  (=  0)  or  (^mamic  (=  1).  The  next  three 
parameters,  NTIME,  NSTEP  and  NOZERO,  are  parameters  for  the  dynamic  analysis 
case  and  are  describes  under  that  sectioa 

2.  Shape  of  Plate  (FEM  Mesh) 

This  program  has  two  methods  for  generating  the  plate  mesh.  Rectangular 
plate  meshes  can  be  generated  by  entering  the  number  of  divisions  between  nodal 
points  along  the  x  and  y-axes,  NX  and  NY,  and  the  lengths  of  each  division,  DX(I) 
and  DY(I)  in  the  x  and  y  directions  respectively.  These  divisions  do  not  have  to  have 
uniform  lengths  because  each  length  is  supplied. 

General  plate  shiq)es  require  their  meshes  to  be  entered  by  the  user.  The  user 
must  provide  the  number  of  elements,  NEM,  the  number  of  nodal  points,  NNM,  the 
element  connectivity  of  the  nodal  points,  NOD(Iv0>  coordinates  of  the  nodal 

points,  X(I)  and  Y(I).  The  plate  can  therefore  be  any  user  defined  shape  including 
curved  edges,  but  must  have  a  uniform  thickness.  The  user  is  encouraged  to  use  other 
mesh  generation  programs  to  develop  ^neral  sh^)ed  (non-rectangular)  meshes. 
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3.  Element  Types 

There  are  three  types  of  quadrilateral  isoparametric  elements  in  this  program. 
Each  nodal  point  in  the  element  has  five  degrees  of  freedom  or  generalized 
displacements  (u,  v,  w,  The  three  types  of  elements  are  four-node  linear, 

eight-node  quadratic,  and  nine-node  quadratic  elements  as  shown  in  Figure  1 8. 


4— Node  8— Node  9-Node 

4  3473473 


Figure  18;  Element  types  and  nodal  point  numbers. 

The  four,  eight,  and  nine-noded  elements  produce  element  matrices  of  20x20, 40x40, 
and  45x45,  respectively.  The  order  of  these  matrices  is  defined  by  the  number  of 
degrees  of  freedom  per  node.  lEL  and  NPE  control  the  type  of  element  used  as 
described  in  section  1. 

The  interpolation  functions  used  for  these  elements  are  isoparametric  and 
belong  to  the  Lagrange  family.  This  allows  the  element  sides  to  be  non-straight  (for 
quadratic  elements).  Full-integration  is  used  for  all  stiffness  terms  except  for  those 
terms  involving  traverse  shear  coefficients  (A44,  A4,,  A,,)  in  which  a  reduced- 
integration  scheme  is  used.  The  reduced-integration  is  performed  to  prevent  shear¬ 
locking  effects. 
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4.  Types  of  Materials  and  Stacking  Sequence 

This  {»rogram  was  developed  for  general  laminated  composite  plates  including 
hybrid  composites.  The  user  suf^lies  the  number  of  different  composite  materials, 
NMTL  in  the  laminate,  and  the  material  properties  of  each  material  (I); 

Elastic  moduli  -  E,  -  E1(I),  Ej  -  E2(I) 

Shear  moduli  -  G.^  -  0120),  G,3  -  G130).  G^  -  0230) 

Poisson's  Ratio  -  v,2  -  ANU120) 

Material  density  •  p  •  RHOO) 

Laminated  plates  have  the  following  characteristics.  The  laminate  is  formed 
by  stacking  a  number  of  orthotropic  material  layers,  NLAY,  in  a  desired  sequence. 
Each  layer,  L,  can  have  different  principal  material  orientations  defined  in  degrees 
from  the  x-axis,  THETAO.^),  be  made  of  different  materials  defined  by  the  material 
number  I,  MH/L),  and  have  varied  diicknesses,  THO^).  The  {H^ogram  is  also  able  to 
analyze  isotropic  materials,  orthotropic  materials  and  hybrids  made  of  these  materials. 

5.  Displacement  on  Boundaries  and  Loading  Conditions 

This  program  accommodates  general  boundary  conditions  by  allowing  the  user 
to  specify  the  number  of  displacement  conditions,  NBDY,  the  location  and  direction 
of  the  displacement  (u,  v,  w,  \|/y),  IBDY(I),  and  the  corresponding  value  of  the 
generalized  displacements,  VBDY(I),  for  any  nodal  point 

There  are  two  ways  to  apply  loads  to  the  plate.  A  uniformly  distributed 
transverse  load  can  be  applied  by  specifying  the  magnitude  and  direction,  PO 
(positive-upward,  negative-downward),  of  the  jnessure.  Other  loads  are  ^rplied  by 
specifying  the  number  of  specified  forces,  NSBF,  the  nodal  point  location  and 
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^neralized  direction,  IBSF(I),  and  magnitude  of  the  load  at  the  nodal  points,  VBSF(I), 
in  the  form  (Nx,  Ny,  Qz,  Mx,  My)  corresponding  to  the  generalized  displacements. 

By  default,  all  generalized  displacements  are  assumed  to  be  fiee  to  move  and 
all  forces  are  assumed  to  be  zero.  The  displacement  and  force  vectors  each  have  five 
times  the  number  of  nodes  entries.  The  ordering  of  the  displacement  and  force 
vectors  are  shown  below: 

Displacement  Vector  -  {{u,v,w,\};,  ,\}>y  }“**  ,{u,v,  w.v*;,  ,v|;y  }“*“ ,...} 

Force  Vector  -  {{n,  ,Ny  ,Q,M,  ,M,  ,{n,  .N,  ,Q.M,  ,M,  (B-1 ) 

To  change  a  boundary  condition  from  the  default  conditions,  the  user  must 
supply  the  position  of  the  displacement  or  force  in  its  respective  vector  and  the  value 
of  that  conditioa  For  example,  to  specify  a  displacement  of  10  in  the  y-direction  on 
node  3,  the  user  would  give  12  for  its  position  in  the  vector  and  10.0  for  its  value. 

One  of  each  of  the  following  boundary  condition  pairs  must  be  specified  along 
the  plate  edges : 

N„oru.,  N„oru,,  Q,  orw,  M„orM;.,  M„or\|;,  (B-2) 

The  following  list  describes  the  applicable  force  and  displacement  values  for 
some  commonly  used  boundary  conditions.  The  coordinates  n  and  s  refer  to  the 
outward  normal  and  the  in-plane  tangential  axes  along  the  plate  edges  respectively. 
For  example  on  an  an  x-edge  (x  =  a),  n  corresponds  to  the  x-axis  and  s  corresponds  to 
the  y-axis. 
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1.  Simply-Supported  Edge 

N.  =  0,  =  0,  w  =  0.  M„  =  0.  V,  =  0  (B-3) 

2.  Hinged  Edge  -  Free  in  the  normal  direction 

N„  =  0,  u,  =  0,  w  =  0,  M,  =  0,  v,  =  0  (B-4) 

3.  Hinged  Edge  -  Free  in  the  tangential  direction 

u,  =  0,  N„  =  0,  w  =  0,  M,  =  0.  v,  =  0  (B-5) 

4.  Clamped  Edge 

u.  =  0,  u.  =  0,  w  =  0,  M/.  =  0,  V,  =  0  (B-6) 

5.  Free  Edge 

N.=0,  N„=0.  Q.  =  0.  M.=0,  M„=0  (B-7) 

6.  Line  of  Symmetry  (for  symmetrical  finite  element  problems) 

u„=0,  N^=0,  Q.=0.  M„=0  (B-8) 


6.  Dynamic  Analysis 

This  program  can  analyze  many  different  dynamic  cases.  For  the  general 
dynamic  case  (ITEM=1),  the  user  must  supply  the  number  of  time  steps,  NTIME,  the 
time  step  size,  DT,  alpha  from  Nevvmark's  method,  AUA,  and  tine  time  step  at  v^ch 
the  transverse  pressure  is  removed,  NSTP.  Additionally,  initial  generalized 
displacements  and  generalized  velocities  can  be  specified  to  add  to  the  range  of 
dynamic  cases  covered  by  the  program.  In  this  case,  all  the  values  of  the  initial 
displacements,  GF0(I),  and  velocities,  GFl(O)  of  each  nodal  point  must  be  entered  in 
the  following  order: 

Displacement  Vector  -  {{u.  v,w,v,,\|;y}“^*,{u,v,w,\}/,,\|/y}“^,... 

Velocity  Vector  -  {{<.  v>,V,,Vy}"^',{u.v,w,\i/,,\i;y}“*“,...} 


(B-9) 


89 


Computer  Program  Structure  (CXKMDPLATE) 

Figure  19  shows  the  general  program  structure  of  COMPLATE  with  all 
subroutine  calls.  The  dashed  boxes  denote  optional  features  of  the  program  depending 
on  control  parameters.  This  is  provided  as  a  reference  for  understanding  the  source 
code. 


MAIN  PROGRAM 


SUBROUTINES 


Open  Data  File 
1— 


Read  Mesh 


X 


Mesh  Generation  (MESH)  | 


Read  Other  Data 


PRE-PROCESSOR 


Compute  Material  Propeties 
I 


Material  Props  (MATPROP) 


Print  Mesh/Material  Irrfb 


I 


Compute  Half-Band  Width 


Compute  Element  Matrices 

. t  . . 


LJ  Element  Stiffness  (STIFF)  LJ  Shape  Functions  (SHAPE) 


Assemble  Global  Matrix 


Add  Force  Conditions 


Add  Displacement  Cond's  |_{Boundary  Conditions(BNDY| 


Solve  Global  Equation  LJ  Equation  Solver  (SOLVE) 


Dynamic  Calculations 


Print  Displacements 


PROCESSOR 


Calculate  Stresses/Strains  u 


Stressea/Strains  (STRESS)  Shape  Functions  (SHAPE) 


Formats 


POST-PROCESSOR 


Figure  19;  Computer  program  structure. 


Data  File  Format 
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The  following  list  describes  the  input  format  for  COMPLATE.  The  data  is 
read  in  using  "free-form"  formal  This  means  the  data  may  be  separated  by  spaces  or 
commas  and  the  only  important  aspect  is  the  order  of  variables  in  the  data  file.  The 
following  list  provides  a  suggested  format  for  ease  in  visually  reading  and  adjusting 
the  input  file.  Descriptions  of  all  variables  used  in  the  program  can  be  found  at  the 
top  of  the  program  source  code  or  at  the  top  of  each  subroutine.  DATA  L1NE(S) 
outlines  how  the  input  is  separated  into  sections  or  lines.  Some  DATA  LINES  require 
more  than  one  actual  input  line  depending  on  the  number  of  variables  to  be  read  in. 
The  information  under  the  TYPE  column  describes  the  data  type  of  the  variable;  Ann 
.  character  string  of  nn  characters,  I  -  integer  (no  decimal  point);  F  -  fixed-point 
(decimal  point  required).  An  example  input  file  is  given  following  this  list  for  clarity. 


VARIABLE 

TYPE 

VARIABLE  DESCRIPTION  AND  NOTES 

DATA  LINE  1 

OUTPUT  HEADING  -  (80  characters) 

TITLE 

A80 

Title  for  ouq)ut  file  -  description  of  problem 

DATA  LINE  2 

PROGRAM  PARAMETERS 

TFT. 

I 

Element  type  (1  =  four  node,  2  =  eight  or  nine  node) 

NPE 

I 

Nodes  per  element  (4  if  IEL=1, 8  or  9  if  IEL=2) 

IMESH 

I 

Indicator  for  rectangular  mesh  generation 
(0  -  all  element  information  is  read  in 

1  -  rectangular  mesh  is  generated) 

NPRNT 

I 

Indicator  for  printing  the  element  stiffness  matrices 
and  force  vectors  (0  -  no  printing,  1  -  printing) 

ITEM 

I 

Indicator  for  transient  analysis  (0  -  static,  1  -  transient) 

NTIME 

I 

Total  number  of  time  steps  (0  for  ITEM  =  0) 

NSTP 

I 

Time  step  number  at  which  loading  is  removed 

NOZERO 

I 

Indicator  for  initial  transient  load  conditions 
(0  -  zero  initial  load,  1  -  non-zero  initial  load) 

SKIP  LINES  3, 4,  AND  5  IF  IMESH  =  1  (MESH  GENERATED) 


DATA  LINE  3 

NEM  I 

NNM  I 

DATALINE(SU 

NOD(I,J)  F 

DATALINETS^S 


READ  MESH  PARAMETERS 

Number  of  elements  in  the  mesh 
Number  of  nodes  in  the  mesh 

MESH  CONNECTIVITY  -  NEM  lines,  NPE  per  line 

Connectivity  of  I-th  element  (J  =  l.NPE) 

NODAL  COORDINATES 


X(I),  Y(I)  F  Global  coordinates  of  I-th  node 

SKIP  LINES  6. 7.  AND  8  IF  IMESH  =  0  (GENERAL  MESH)  •** 


DATA  LINE  6 


MESH  GENERATION  PARAMETERS 


NX  I 

NY  I 

DATALINE(S^7 

DX(I)  F 

DATALINE(S)8 

DY(I)  F 

DATA  LINE  9 

NMTL  I 


Number  of  element  subdivisions  along  the  x-axis 
Number  of  element  subdivisions  along  the  y-axis 

X-DIVISIONS  -  IEL*NX  Entries 

Distance  between  two  nodes  along  the  x-axis 

Y-DIVISIONS  -  IEL*NY  Entries 

Distance  between  two  nodes  along  the  y-axis 

MATERIALS 

Number  of  different  materials  in  the  laminate 


DATALINEfSUO 


E1(I)  F 

E2(I)  F 

G12(I)  F 

G13a)  F 

G23(I)  F 

ANU12(I)  F 

RHO(I)  F 


MATERIAL  I  PROPERTIES  -  NMTL  lines 
Modulus  along  fiber  direction  (1 -direction) 

Modulus  transverse  to  fiber  direction  (2-direction) 
In-plane  shear  modulus  oriented  along  fiber  direction 
Shear  modulus  with  respect  to  1-3  plane 
Shear  modulus  with  respect  to  2-3  pl»ne 
In-plane  Poisson's  ratio  (1-2  plane) 

Material  Density 


DATA  LINE  11 

NLAY  I 

DATALINEIS)  12 

MTL(I)  I 

THETA(I)  F 

TH(I)  F 

DATA  LINE  13 

PO  F 

DATA  LINE  14 

NBDY  I 

DATALINE(SH5 

IBDY(I)  I 

DATALINEISl  16 
VBDY(I)  F 

DATA  LINE  17 
NBSF  I 
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NUMBER  OF  LAMINA 

Number  of  lamina  (layers)  in  the  laminate 

LAMINA  PROPERTIES  -  NLAY  lines 

Material  number  of  I-th  lamina 

Fiber  orientation  angle  of  1-th  lamina  (in  degrees) 

Thickness  of  I-th  lamina 

UNIFORM  TRANSVERSE  LOADING 

Intensity  of  uniformly  distributed  transverse  load 

SPECIFIED  DISPLACEMENTS  (cannot  be  zero) 

Number  of  specified  generalized  displacements 

DISPLACEMENTS  -  NBDY  entries 

Location/direction  of  specified  displacement  I 
Order  -  by  node  number  and  (u,  v,  w,  Vy) 

DISPLACEMENT  VALUES  -  NBDY  entries 

Value  of  displacement  corresponding  to  IBDY(I) 

SPECIFIED  FORCES 

Number  of  specified  generalized  forces 


SKIP  LINES  18  AND  19  IF  NBSF  =  0  (NO  FORCE  CONDITIONS  OTHER 

THAN  TRANSVERSE  PRESSURE)  *** 


FORCES  -  NBSF  entries 


IBSF(I)  I  Location/direction  of  specified  generalized  forces 

Order  -  by  node  number  and  (N^,  Ny,  Q2,  M^,  My) 


VBSF(I) 


F 


FORCE  VALUES  -  NBSF  entries,  8  per  line 
Value  of  specified  force  corresponding  to  1BSF(I) 


**♦  SKIP  LINES  20, 21,  AND  22  IF  ITEM  =  0  (STATIC  ANALYSIS) 
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Sample  Input  Data  File 


STATIC 

BENDING  OF 

A  SIMPLY 

SUPPORTED 

PLATE  UNDER  UNIFORM 

LOAD  (2X2Q9  MESH) 

2 

9 

1 

0  0 

1 

0 

0 

2 

2 

1.25 

1.25 

1.25 

1.25 

1.25 

1 

25.  E6 

*5 

1.25 

1.25 

1.25 

1.E6 

0.5E6 

0.5E6 

0.2E6 

0.25 

0.3 

0 

1 

0.0 

0.033333 

1 

90.0 

0.033334 

1 

0.0 

0.033333 

1.0 

37 

1 

2 

4 

5  7 

10 

12 

15  17 

20  22 

23  25 

26 

29 

48 

50 

51 

54  73 

75 

76 

79  98 

100  101 

103  104 

108 

109 

113 

114 

118 

119  123 

124 

125 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0 


Sample  Output  File 


The  following  output  conesponds  to  the  previous  input  Hie. 


STATIC  BENDING  OF  A  SIMPLY  SUPPORTED  PLATE  UNDER  UNIFORM  LOAD  (2X209  MESH) 

ELEMENT  TYPECl-LINEAR. 2-QUADRATIC)  -  2  NODES  PER  ELEMENT-  9 
ACTUAL  NUMBER  OF  ELEMENTS  IN  THE  MESH-  4 
NUMBER  OF  NODES  IN  THE  MESH  -  25 
DEGREES  OF  FREEDOM  -  5 

MATERIAL  1  PROPERTIES; 

MODULUS. El-  0.25000E+08 
MODULUS. E2-  O.lOOOOE+07 
SHEAR  MODULUS. G12-  0.50000E+06 
SHEAR  MODULUS. G13-  0.50000E+06 
SHEAR  MODULUS. G23-  0.20000E+06 
POISSONS  RATIO. NU12-  0.25000E+00 
MATERIAL  DENSITY. RHO-  0.30000E+00 


LAMINATE  STACKING  SEOUENCE 


LAYER  MTL  # 

1  1 

2  1 

3  1 


THETA  THICKNESS 
0.00000  0.33333E-01 
90.00000  0.33334E-01 
0.00000  0.33333E-01 


TOTAL  THICKNESS  -  O.lOOOOE+00 
LAMINATE  PLATE  PROPERTIES 


A  MATRIX  TERMS 
0.17042E+07 
0.25063E+05 
O.OOOOOE+00 

0.25063E+05 

0.90227E+06 

O.OOOOOE+00 

O.OOOODE+00 

O.OOOOOE+00 

0.50000E+05 

B  MATRIX  TERMS 
O.OOOOOE+00 
O.OOOOOE+00 
O.OOOOOE+00 

O.OOOOOE+00 

O.OOOOOE+00 

O.OOOOOE+00 

O.OOOOOE+00 

O.OOOOOE+00 

O.OOOOOE+00 

D  MATRIX  TERMS 
0.20143E+04 
0.20886E+02 
O.OOOOOE+00 

0. 20886 E+02 
0.16781E+03 
O.OOOOOE+00 

O.OOOOOE+DO 

O.OOODOE+00 

0.41667E+02 

SHEAR  TERMS:  A44. 
0.25000E+05 

A45.  ASS 
O.OOOOOE+00 

0.33333E+05 

INERTIAL  TERMS  11 
0.30000E-01 

.  12.  13 
O.OOOOOE+00 

0.25000E-04 

NUMBER  OF  SPECIFIED 

DISPLACEMENTS- 
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SPECIFEO  DISPLACEMENTS  AND  THEIR  VALUES  FOLLOW: 


1 

O.OOOOOE+00 

2 

O.OOOOOE+00 

4 

O.OOOOOE+OO 

6 

O.OOOOOE+00 

7 

O.OOOOOE+00 

10 

U.OOOOOE+00 

12 

O.OOOOOE+00 

15 

O.OOOOOE+00 

17 

O.OOOOOE+OO 

20 

O.OOOOOE+00 

22 

O.OOOOOE+00 

23 

O.OOOOOE+OO 

25 

O.OOOOOE+00 

26 

O.OOOOOE+OO 

29 

O.OOOOOE+OO 

48 

O.OOOOOE+00 

50 

O.OOOOOE+00 

51 

O.OOOOOE+OO 

54 

O.OOOOOE+00 

73 

O.OOOOOE+OO 

75 

O.OOOOOE+OO 

76 

O.OOOOOE+00 

79 

O.OOOOOE+OO 

98 

O.OOOOOE+OO 

100 

O.OOOOOE+00 

101 

O.OOOOOE+OO 

103 

O.OOOOOE+OO 

104 

O.OOOOOE+00 

108 

O.OOOOOE+OO 

109 

O.OOOOOE+OO 

113 

O.OOOOOE+00 

114 

O.OOOOOE+OO 

118 

O.OOOOOE+OO 

119 

125 

o  o 

o  o 
o  o 
o  o 
o  o 
o  o 
m  m 
+  + 
O  O 

o  o 

123 

O.OOOOOE+OO 

124 

O.OOOOOE+OO 

UNIFORMLY  DISTRIBUTED  LOAD.  PO  -  O.lOOOOE+01 
NUMBER  OF  SPECIFIED  FORCES-  0 

SPECIFIED  FORCE  DEGREES  OF  FREEDOM  AND  THEIR  VALUES  FOLLOW: 


BOOLEAN  (CONNECTIVITY)  MATRIX-NOD{ I .J ) 


1  1  3 

13  11  2 

8 

12  6  7 

2  3  5 

15  13  4 

10 

14  8  9 

3  11  13 

23  21  12 

18 

22  16  17 

4  13  15 

25  23  14 

20 

24  18  19 

COORDINATES  OF  THE 

GLOBAL  NODES: 

1 

O.OOOOOE+OO 

O.OOOOOE+OO 

2 

0.12500E+01 

O.OOOOOE+OO 

3 

0.25000E+01 

O.OOOOOE+OO 

4 

0.37500E+01 

O.OOOOOE+OO 

5 

0.50000E+01 

O.OOOOOE+OO 

6 

O.OOOOOE+OO 

0.12500E+01 

7 

0.12500E+01 

0.12500E+01 

8 

0.25000E+01 

0.12600E+01 

9 

0.37500E+01 

0.12500E+01 

10 

0.50000E+01 

0.12500E+01 

11 

O.OOOOOE+OO 

0.25000E+01 

12 

0.12500E+01 

0.25000E+01 

13 

0.25000E+01 

0.25000E+01 

14 

0.37500E+01 

0.25000E+01 

15 

0.50000E+01 

0.25000E+01 

16 

O.OOOOOE+OO 

0.37500E+01 

17 

0.12500E+01 

0.37500E+01 

18 

0.25000E+01 

0.37500E+01 

19 

0.37500E+01 

0.37500E+01 

20 

0.50000E+01 

0.37500E+01 

21 

O.OOOOOE+OO 

0.50000E+01 

22 

0.12500E+01 

0.50000E+01 

23 

0.25000E+01 

0.50000E+01 

24 

0.37500E+01 

0.50000E+01 

25 

0.50000E+01 

0.50000E+01 

HALF 

BAND  WIDTH  OF 

GLOBAL  STIFFNESS  MATRIX  -  65 

GENERALIZED  DISPLACEMENTS  (U.V.W.SX.SY)  PER  NODE 


1 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.67040E-01 

O.OOOOOE+OO 

O.OOOOOE+OO 

2 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.6201ZE-01 

0.78426E-02 

O.OOOOOE+OO 

3 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.47771E-01 

0.14677E-01 

O.OOOOOE+OO 

4 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.25989E-01 

0.19509E-01 

O.OOOOOE+OO 

5 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.21343E-01 

O.OOOOOE+OO 

6 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.63455E-01 

O.OOOOOE+OO 

0.59374E-02 

7 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.58703E-01 

0.74048E-02 

0. 54801 E- 02 

8 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.45230E-01 

0.13868E-01 

0.41977E-02 

9 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.24611E-01 

0.18440E-01 

0.22724E-02 

10 

O.OOOOOE+00 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.20175E-01 

O.OOOOOE+OO 

11 

O.OOOOOE+00 

O.OOOOOE+OO 

0.51706E-01 

O.OOOOOE+OO 

0.13060E-01 

12 

O.OOOOOE+00 

O.OOOOOE+OO 

0.47869E-01 

0.59961E-02 

0.12039E-01 

13 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.36951E-01 

0.11279E-01 

0.91828E-02 

14 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.20143E-01 

0.15090E-01 

0.49446E-02 

15 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.16560E-01 

O.OOOOOE+OO 

16 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.29794E-01 

O.OOOOOE+OO 

0.21322E-01 

17 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.27620E-01 

0.34017E-02 

0.19733E-01 

18 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.21398E-01 

0.64412E-02 

0.15228E-01 

19 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.11718E-01 

0.87228E-02 

0.82896E-02 

20 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.96390E-02 

O.OOOOOE+OO 

21 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.25620E-01 

22 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.23769E-01 

23 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.18493E-01 

24 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

0.10180E-01 

25 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 

LAMINATE  STRAINS  AND  STRESSES  AT  GAUSS  POINTS 
(X-COORD.Y-COORO) 

LOG  EPSILONX  EPSILONY  GAMMAXY  KAPPAX  KAPPAY  XAPPAXY 
LAY  Z- COORD  SIGMAX  SIGMAY  SIGMAXY  SIGMAYZ  SI6MAXZ 

(  0.5283E+00.  0.5283E+00) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.6283E-02  D.4604E-02  •0.2520E-03 

BOT  -0.3141E-03  -0.2302E-03  0.1260E-04 

TOP  0.3141E-03  0.2302E-03  -0.1260E-04 

1  -0.5000E-01  -0.7931E+04  -0.3095E+03  0.63006*01  0.41256*00  -0.8524E+01 

1  -0.1667E-01  -0.2644E+04  -0.1032E+03  0.2100E+01  0.4125E+00  ■0.8524E+0I 

2  -0.1667E-01  -0.1242E+03  ■0.1949E+04  0.2100E+01  0.1031E+01  -0.3410E+01 

2  0.1667E-01  0.1242E+03  0.1949E+04  -0.2100E+01  0.1031E+01  -0.3410E+01 

3  0.1667E-01  0.2644E+04  0.1032E+03  -0.2100E+01  0.4125E+00  -0.8524E+01 

3  0.5000E-01  0.7931E+04  0.3095E+03  -0.6300E+01  0.4125E+00  •0.8524E+01 

(  0.5283E+00.  0.1972E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.5419E-02  0.5678E-02  -0.1055E-02 

BOT  -0.2710E-03  -0.2839E-03  0.5277E-04 

TOP  0.2710E-03  0.2839E-03  -0.5277E-04 

1  -0.5000E-01  -0.6862E+04  -0.3525E+03  0.2639E+02  0.4592E+00  -0.7366E+01 

1  -0.1667E-01  -0.2287E+04  -0.1175E+03  0.8795E+01  0.4592E+00  -0.7366E+01 

2  -0.1667E-01  -0.1143E+03  -0.2394E+04  0.8795E+01  0.1148E+01  -0.2946E+01 

2  0.1667E-01  0.1143E+03  0.2394E+04  -0.8795E+01  0.1148E+01  -0.2946E+01 

3  0.1667E-01  0.22876*04  0.1175E+03  -0.8795E+01  0.4592E+00  -0.7366E+01 

3  0.5000E-01  0.6862E+04  0.3525E+03  -0.2639E+02  0.4592E+00  -0.7366E+01 

(  0.1972E+01.  0.5283E+00) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.5359E-02  0.3815E-02  -0.8633E-03 

BOT  -0.2680E-03  -0.1907E-03  0.4316E-04 

TOP  0.2680E-03  0.1907E-03  -0.4316E-04 

1  -0.5000E-01  -0.6763E+04  -0.2584E+03  0.2158E+02  0.3032E+00  -0.3156E+02 

1  -0.1667E-01  -0.2255E+04  -0.8612E+02  0.7194E+01  0.3032E+00  -0.3156E+02 

2  -0.1667E-01  -0.1055E+03  -0.1616E+04  0.7194E+01  0.7579E+00  -0.1263E+02 

2  0.1667E-01  0.1055E+03  0.1616E+04  -0.7194E+01  0.7579E+00  -0.1263E+02 

3  0.1667E-01  0.2255E+04  0.8612E+02  -0.7194E+01  0.3032E+00  -0.3156E+02 

3  0.5000E-01  0.6763E+04  0.2584E+03  -0.2168E+02  0.3032E+00  -0.3156E+02 
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(  0.1972E+01.  0.1972E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.4653E-02  0.4676E-02  -0.3595E-02 

BOT  -0.2327E-03  -0.2338E-03  0.1797E-03 

TOP  0.2327E-03  0.2338E-03  -0.1797E-03 

1  -0.5000E-01  -0.5890E+04  -0.2927E+03  0.8987E+02  0.4438E+00  -0.2798E+02 

1  -0.1667E-01  -0.1963E+04  -0.9758E+02  0.2996E+02  0.4438E+00  -0.2798E+02 

2  -0.1667E-01  -0.9728E+02  -0.1973E+04  0.2996E+')2  O.lllOE+Ol  -0.1119E+02 

2  0.1667E-01  0.9728E+02  0.1973E+04  -0.2996E+02  O.lllOE+01  -0.1119E+02 

3  0.1667E-01  0.1963E+04  0.9758E+02  -0.2996E+02  0.4438E+00  -0.2798E+02 

3  0.5000E-01  0.5890E+04  0.2927E+03  -0.8987E+02  0.4438E+00  ■0.2798E+02 

(  0.3028E+01,  0.5283E+00) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.4011E-02  0.2702E-02  -0.1229E-02 

BOT  -0.2005E-03  -0.1351E-03  0.6144E-04 

TOP  0.2005E-03  0.1351E-03  -0.6144E-04 

1  -0.5000E-01  -0.5060E+04  -0.1857E+03  0.3072E+02  0.1910E+00  -0.4811F+02 

1  -0.1667E-01  -0.1687E+04  -0.6191E+02  0.1024E+02  0.1910E+00  -0.4811E+02 

2  -0.1667E-01  -0.7830E+02  -0.1146E+04  0.1024E+02  0.4775E+00  -0.1924E+02 

2  0.1667E-01  0.7830E+02  0.1  46E+04  -0.1024E+02  0.4775E+00  -0.1924E+02 

3  0.1667E-01  0.1687E+04  0.6191E+02  -0.1024E+02  0.1910E+00  -0.4811E+02 

3  0.5000E-01  0.5060E+04  0.1857E+03  -0.3072E+02  0.1910E+00  -0.4811E+02 

(  0.3028E+01.  0.1972E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.3516E-02  0.3280E-02  -0.5039E-02 

BOT  -0.1758E-03  -0.1640E-03  0.2520E-03 

TOP  0.1758E-03  0.1640E-03  -0.2520E-03 

1  -0.5000E-01  -0.4447E+04  -0.2085E+03  0.1260E+03  0.1813E+00  -0.4337E+02 

1  -0.1667E-01  -0.1482E+04  -0.6949E+02  0.4200E+02  0.1813E+00  -0.4337E+02 

2  -0,1667E-01  -0.72442+02  -0.1385E+04  0.4200E+02  0.4533E+00  -0.1735E+02 

2  0.1667E-01  0.7244E+02  0.1385E+04  -0.4200E+02  0.4533E+00  -0.1735E+02 

3  0.1667E-01  0.1482E+04  0.6949E+02  -0.4200E+02  0.1813E+00  -0.43372+02 

3  0.5000E-01  0.4447E+04  0.2085E+03  -0.1260E+03  0.1813E+00  -0.4337E+02 

(  0.4472E+01.  0.5283E+00) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.1266E-02  0.7917E-03  -0.1505E-02 

BOT  -0.6329E-04  -0.3958E-04  0.7524E-04 

TOP  0.6329E-04  0.3958E-04  -0.7524E-04 

1  -0.5000E-01  -0.1596E+04  -0.5564E+02  0.3762E+02  0.1455E+00  -0.7036E+02 

1  -0.1667E-01  -0.5321E+03  -0.1852E+02  0.1254E+02  0.1465E+00  -0.7036E+02 

2  -0.1667E-01  -0.2446E+02  -0.3360E+03  0.1254E+02  0.3639E+00  -0.2814E+02 

2  0.1667E-01  0.2446E+02  0.3360E+03  -0.1254E+02  0.3639E+00  -0.2814E+02 

3  0.1667E-01  0.5321E+03  0.1852E+02  -0.1254E+02  0.1455E+00  -0.7036E+02 

3  0.5000E-01  0.1596E+04  0.5554E+02  -0.3762E+02  0.1455E+00  -0.7036E+02 

(  0.4472E+01,  0.1972E+01) 

MID  O.OOOOE+00  0.0000E+.  O.OOOOE+00  0.1121E-D2  0.9492E-03  -D.6048E-02 

BOT  -0.6606E-04  -0.4746E-04  0.3024E-03 

TOP  0.5606E-04  0.4746E-04  -0.3024E-03 

1  -0.5000E-01  -0.1417E+04  -0.6163E+02  0.1612E+03  0.2988E-01  -0.6418E+02 

1  -0.1667E-01  -0.4723E+03  -0.2054E+02  0.5040E+02  0.2988E-01  -0.6418E+02 

2  -0.1667E-01  -0.2270E+02  -0.4012E+03  0.5040E+02  0.7469E-01  -0.2567E+02 

2  0.1667E-01  0.2270E+02  0.4012E+03  -0.5040E+02  0.7469E-01  -0.2567E+02 

3  0.1667E-01  0.4723E+03  0.2054E+02  -0.5040E+02  0.2988E-01  -0.6418E+02 

3  0.5000E-01  0.1417E+04  0.6163E+02  -0.1512E+03  0.2988E-01  -0.6418E+02 

(  0.5283E+00,  0.3028E+01) 

MID  O.OOOOE+00  O.OOOOE+DO  O.OOOOE+00  0.4D34E-02  D.6765E-D2  -0.1782E-02 

BOT  -0.2017E-03  -0.3382E-03  0.8912t-04 
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TOP  0.2017E-03  0.3382E-03  -0.8912E-04 

1  -0.5000E-01  -0.5139E+04  -0.3896E-t-03  0.4456E+02  -0.1323E+01  -0.4661E+01 

1  -0.1667E-01  -0.1713E+04  -0. 1299E+03  0.1485E+02  -0.1323E+01  -0.466^+01 

2  -0.1667E-01  -0.9565E+02  -0.2843E+04  0.1485E-t-02  -0.3306E+01  -0.1865E+01 

2  0.1667E-01  0.9565E+02  0.2843E+04  -0.1485E+02  -0.3306E+01  -0.1865E+01 

3  0.1667E-01  0.1713E+04  0.1299E+03  -0.1485E+02  -0.1323E+01  -0.4661E+01 

3  0.5000E-01  0.5139E+04  0.3896E+03  -0.4456E+02  -0.1323E+01  -0.4661E+01 

(  0.5283E+00.  0.4472E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.1239E-02  0.3158E-02  -0.2425E-02 

BOT  -0.6193E-04  -0.1579E-03  0.1213E-03 

TOP  0.6193E-04  0.1579E-03  ■0.1213E03 

1  -C.5000E-01  -0.1592E+04  -0.1738E+03  0.6063E+02  -0.8194E+01  -0.1608E+01 

1  -0.1667E-01  -0.5306E+03  -0.5794E+02  0.2021E+02  -0.8194E+01  -0.1608E+01 

2  -0.1667E-01  -0.3389E+02  -0.1324E+04  0.2021E+02  -0.2049E+02  -0.6432E+00 

2  0.1667E-01  0.3389E+02  0.1324E+04  -0.2021E+02  ■0.2049E+02  ■0.6432E+00 

3  0.1667E-01  0.5306E+03  0.5794E+02  -0.2021E+02  -0.8194E+01  -0.1608E+01 

3  0.5000E-01  0.1592E+04  0.1738E+03  -0.6063E+02  -0.8194E+01  -0.1608E+01 

(  0.1972E+01.  0.3028E+01) 

MID  O.OOOOE+00  0. 0000^+00  O.OOOOE+00  0.3510E-02  0.5699E-02  -0.6142E-02 

BOT  -0.1755E-03  -0.2849E-03  0.3071E-03 

TOP  0.1755E-03  0.2849E-03  -0.3071E-03 

1  -0.5000E-01  -0.4470E+04  -0.3297E+03  0.1535E+03  -0.9461E+00  -0.1815E+02 

1  -0.1667E-01  -0.1490E+04  -0.1099E+03  0.5118E+02  -0.9461E+00  -0.1815E+02 

2  -0.1667E-01  -0.8246E+02  -0.2395E+04  0.5118E+02  -0.2365E+01  -0.7258E+01 

2  0.1667E-01  0.8246E+02  0,2395E+04  -0.5118E+02  ■0.2365E+01  -0.7268E+01 

3  0.1667E-01  0.1490E+04  0.1099E+03  -0.6118E+02  ■0.9461E+00  -0.1815E+02 

3  0.5000E-01  0.4470E+04  0.3297E+03  -0.15352+03  -0.9461E+00  -0.1815E+02 

(  0.1972E+01.  0.4472E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.1096E-02  0.2723E-D2  -0.8488E-02 

BOT  -0.5479E-04  -0.1361E-03  0.4244E-03 

TOP  0.5479E-04  0.1361E-03  -0.4244E-03 

1  -0.5000E-01  -0.1407E+04  -0.1502E+03  0.2122E+03  -0.7229E+01  -0.6760E+01 

1  -0.1667E-01  -0.4691E+03  -0.5007E+02  0.7074E+02  -0.7229E+01  -0.676QE+01 

2  -0.1667E-01  -0.2968E+02  -0.1142E+04  0.7074E+02  -0.1807E+02  -0.2704E+C1 

2  0.1667E-01  0.2968E+02  0.1142E+04  -0.7074E+02  -0.1807E+02  -0.2704E+01 

3  0.1667E-01  0.4691E+03  0.5007E+02  -0.7074E+02  -0.7229E+01  -0.6760E+01 

3  0.5000E-01  0.1407E+04  0.1502E+03  -0.2122E+03  -0.7229E+01  -0.6760E+01 

(  0.3028E+01.  0.3028E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.2727E-02  0.4125E-02  -0.8732E-02 

BOT  -0.1364E-03  -0.2062E-03  0.4366E-03 

TOP  0.1364E-03  0.2062E-03  -0.4366E-03 

1  -0.5000E-01  -0.3470E+04  -0.2409E+03  0.2183E+03  -C.3552E+00  -0.2968E+02 

1  -0.1667E-01  -0.1157E+04  -0.8031E+02  0.7277E+02  -0.3552E+00  -0.2968E+02 

2  -0.1667E-01  -0.6280E+02  -0.1734E+04  0.7277E+02  -0.8881E+00  -0.1187E+02 

2  0.1667E-01  0.6280E+02  0.1734E+04  -0.7277E+02  -0.8881E+00  -0.1187E+02 

3  0.1667E-01  0.n57E+04  0.8031E+02  -0.7277E+02  -0.3552E+00  -0.2968E+02 

3  0.5000E-01  0.3470E+04  0.2409E+03  -0.2183E+03  -0.3552E+00  -0.2968E+02 

(  0.3028E+01,  0.4472E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.8835E-03  0.2059E-02  -0.1239E-01 

BOT  -0.4418E-04  -0.1029E-03  0.6196E-03 

TOP  0.4418E-04  0.1029E-03  -0.6196E-03 

1  -0.5000E-01  -0.1133E+04  -0.1143E+03  0.3098E+03  -0.5192E+01  -0.1112E+02 

1  -0.1667E-01  -0.3777E+03  -0.3809E+02  0.1033E+03  -0.5192E+01  -0.1112E+02 


2  -0.1667E-01  -0.2336E+02  -0.8637E+03  0.1033E+03  -0.1298E+02  -0.4450E+01 

2  0.1667E-01  0.2336E+02  0.8637E+03  -0.1033E+03  -0.1298E+02  -0.4450E+01 

3  0.1667E-01  0.3777E+03  0.3809E+02  -0.1033E+03  -0.5192E+01  ■0.1112E+02 

3  0.5000E-01  0.1133E+04  0.1143E+03  -0.3098E+03  -0.5192E+01  -0.1112E+02 

(  0.4472E+01.  0.3028E+01) 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.9018E-03  0.1233E-02  -0.1057E-01 

BOT  -0.4509E-04  -0.6165E-04  0.5287E-03 

TOP  0.4509E-04  0.6165E-04  -0.5287E-03 

1  -0.5000E-01  -0.1146E+04  -0.7310E+02  0.2643E+03  0.7850E-01  -0.4947E+02 

1  -0.1667E-01  -0.3819E+03  -0.2437E+02  0.8811E+02  0.7850E-01  -0.4947E+02 

2  -0.1667E-01  -0.2022E+02  -0.5188E+03  0.8811E+02  0.1963E+00  -0.1979E+02 

2  0.1667E-01  0.2022E+02  0.5188E+03  -0.8811E+02  0.1963E+00  -0.1979E+02 

3  0.1667E-01  0.3819E+03  0.2437E+02  -0.8811E+02  0.7850E-01  -0.4947E+02 

3  0.5000E-01  0.1146E+04  0.7310E+02  -0.2643E+03  0.7850E-01  -0.4947E+02 

(  0.4472E+01,  0.4472E+0n 

MID  O.OOOOE+00  O.OOOOE+00  O.OOOOE+00  0.3065E-03  0.6503E-03  -0.1545E-01 

BOT  -0.1533E-04  -0.3251E-04  0.7723E-03 

TOP  0.1533E-04  0.3251E-04  -0.7723E-03 

1  -0.5000E-01  -0.3923E+03  -0.3644E+02  0.3861E+03  -0.2025E+01  -0.2162E+02 

1  -0.1667E-01  -0.1308E+03  -0.1215E+02  0.1287E+03  -0.2025E+01  -0.2162E+02 

2  -0.1667E-01  -0.7838E+01  -0.2729E+03  0.1287E+03  -0.5063E+01  -0.8647E+01 

2  0.1667E-01  0.7838E+01  0.2729E+03  -0.1287E+03  -0.5063E+01  -0.8647E+01 

3  0.1667E-01  0.1308E+03  0.1215E+02  -0.1287E+03  -0.2025E+01  -0.2162E+02 

3  0.5000E-01  0.3923E+03  0.3644E+02  -0.3861E+03  -0,2025E+01  -0.2162E+02 


APPENDIX  C 


COMPLATE  -  COMPUTER  PROGRAM  SOURCE  CODE 


C  COMPUTER  PROGRAM  COMPLATE 

C  (STATIC  AND  TRANSIENT  ANALYSIS  OF  COMPOSITE  PLATES) 

C 

C  Revision  of  program  PLATE  by  J.N.  Reddy  for  orthotropic  plates 
C  Source:  Reddy,  J.N.  'Introduction  to  the  Finite  Element  Method*. 

C  McGraw-Hill.  New  York  (1984). 

C 

C  Revision  by:  Brett  A.  Pauer.  The  Ohio  State  University 
C 

C  . 

C  . 


C  .  DESCRIPTION  OF  THE  VARIABLES 

C  . 

C  .  A(I,J) . EXTENSIONAL  STIFFNESS  MATRIX  (I.J-1.2.3) 

C  .  A(K.L) . SHEAR  TERMS  OF  STIFFNESS  MATRIX  (K,L-4,5) 

C  .  AO. A1.A2. A3. A4... PARAMETERS  IN  THE  TIME-APPROXIMATION  SCHEME 

C  .  AK . SHEAR  CORRECTION  COEFFICIENT 

C  .  ALFA . PARAMETER  IN  THE  NEWMARK  SCHEME 

C  ,  ANU12(I) _ POISSON’S  RATIO  OF  MATERIAL 

C  .  B(I,J) . COUPLING  STIFFNESS  MATRIX  (I,J-1.2,3) 

C  ,  BETA . PARAMETER  IN  THE  NEWMARK  SCHEME 

C  .  D(I.J) . BENDING  STIFFNESS  MATRIX  (I.J-1.2.3) 

C  .  DX( I ).DY( I). DISTANCE  BETWEEN  NODES  IN  X.Y  DIRECTIONS  FOR  MESH 
C  .  GENERATION 

C  .  DT . TIME  INCREMENT  IN  THE  TRANSIENT  ANALYSIS 

C  .  Eld). E2(I). ELASTIC  MODULI  OF  MATERIAL  I 
C  .  ELP(I) . ELEMENT  FORCE  VECTOR 

C  .  ELXY( I. J)... ELEMENT  NODE  COORDINATES  OF  ELEMENT  NODE  I 
C  .  J-1  FOR  X- COORD.  J-2  FOR  Y- COORD 

C  .  G12( I ).G23( I ).G13( I).. SHEAR  MODULI  OF  MATERIAL  I 

C  .  6F(I) . GLOBAL  FORCE  VECTOR;  SOLUTION  VECTOR  FROM  ’SOLVE’ 

C  .  GFO(I) . SOLUTION  VECTOR  AT  CURRENT  TIME 

C  .  GFKI) . FIRST  TIME  DERIVATIVE  OF  THE  SOLUTION  (VELOCITY) 

C  .  GF2(I) . SECOND  TIME  DERIVATIVE  OF  THE  SOLUTION  (ACCELERATION) 

C  .  GSTIF(N.M).. GLOBAL  STIFFNESS  MATRIX  (IN  BANDED  FORM) 

C  .  H . THICKNESS  OF  THE  PUTE 

C  .  IBDY(I) . LOCATION  AND  DIRECTION  OF  SPECIFIED  GLDBAL 

C  .  DISPLACEMENTS 

C  .  IBSF(I) . LOCATION  AND  DIRECTION  OF  SPECIFIED  NONZERO  GLOBAL 

C  .  FORCES 

C  .  lEL . INDICATOR  FOR  THE  ELEMENT  TYPE: 

C  .  I  EL-1,  4 -NODE  ELEMENT 

C  .  IEL-2.  8-  OR  9-NODE  ELEMENT 


100 
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C  .  IMESH . INDICATOR  FOR  MESH  GENERATION 

(0-READ  IN.  1 -SQUARE  MESH  IS  GENERATED) 

ITEM . INDICATOR  FOR  TRANSIENT  ANALYSIS  (1-YES,  O-NO) 

MTL(L) . MATERIAL  NUMBER  OF  LAYER  L 

NBOY . TOTAL  NUMBER  OF  SPECIFIED  DEGREES  OF  FREEDOM 

NBSF . TOTAL  NUMBER  OF  SPECIFIED  NONZERO  FORCES 

NCMAX . VALUE  OF  THE  COLUMN- DIMENSION  OF  GST  IF 

NDF . NUMBER  OF  DEGREES  OF  FREEDOM  PER  NODE  (U.V.H.SX.SY) 

NEM . NUMBER  OF  ELEMENTS 

NEO . TOTAL  NUMBER  OF  DEGREES  OF  FREEDOM  (NODESxNODAL  DOF) 

NHBW . HALF-BAND  WIDTH  OF  GLOBAL  STIFFNESS  MATRIX 

NLAY . NUMBER  OF  LAYERS  IN  THE  UMINATE 

NMTL . NUMBER  OF  DIFFERENT  MATERIALS  IN  THE  LAMINATE 

NN . NUMBER  OF  DEGREES  OF  FREEDOM  PER  NODE 

(NODES  PER  ELEMENT  x  NODAL  DOF) 

NNM . NUMBER  OF  GLOBAL  NODES 

NOD(I.J) _ ELEMENT  CONNECTIVITY  MATRIX 

NOZERO . INDICATOR  FOR  ZERO{NOZERO-0)  OR  NONZERO! NOZERO-1 ) 

INITIAL  CONDITIONS  FOR  TRANSIENT  ANALYSIS 

NPE . NUMBER  OF  NODES  PER  ELEMENT  (4.  8  OR  9) 

NPRNT . INDICATOR  FOR  PRINTING  ELEMENT  MATRICES  AND  FORCE 

VECTORS  (1-PRINT.  0-00  NOT  PRINT) 

NRMAX . VALUE  OF  THE  ROW-DIMENSION  OF  GSTIF 

NSTP . TIME  STEP  AT  WHICH  THE  LOAD  IS  REMOVED  FROM  THE 

PLATE  (IN  THE  TRANSIENT  ANALYSIS) 

NT . CURRENT  TIME  STEP  NUMBER  IN  THE  TRANSIENT  ANALYSIS 

NTIME . NUMBER  OF  TIME  STEPS  IN  THE  TRANSIENT  ANALYSIS 

NX . NUMBER  OF  DIVISIONS  ALONG  X-AXIS  FOR  MESH  GENERATION 

NY . NUMBER  OF  DIVISIONS  ALONG  Y-AXIS  FOR  MESH  GENERATION 

PO . INTENSITY  OF  APPLIED  TRANSVERSE  UNIFORM  PRESSURE 

QBAR( I. J.L), TRANSFORMED  STRESS-STRAIN  MATRIX  OF  LAYER  L 

RHO(I) . DENSITY  OF  MATERIAL  I 

STIF(N.M)... ELEMENT  STIFFNESS  MATRIX 

T . TIME  VARIABLE  IN  THE  TRANSIENT  ANALYSIS 

TH(L) . THICKNESS  OF  LAYER  L 

THETA(L),... FIBER  DIRECTION  ORIENTATION  OF  LAYER  L 
TITLE . TITLE  FROM  INPUT  DATA  FILE 

VBDY(I) . VALUES  OF  THE  DISPLACEMENTS  CORRESPONDING  TO  IBDY(I) 

VBSF(I) . VALUES  OF  SPECIFIED  FORCES  CORRESPONDING  TO  IBSF(I) 

W0,W1.W2 _ ARRAYS  CORRESPONDING  TO  GF0.GF1.GF2  IN  AN  ELEMENT 

X(I).Y(I)...X  AND  Y  COORDINATES  OF  GLOBAL  NODE  I 


IMPLICIT  REAL*8(A-H.0-Z) 

CHARACTER  DATFILE*20,0UTPTF*20 

DIMENSION  GSTIF(1000.200),GF(500).GF0(500).GF1(500).GF2(500). 

*  VBOY(400).IBOY(400).VBSF(400).IBSF(400).TITLE(20). 

*  E1(20),E2(20),GI2(20),G23(20).G13(20).ANU12(20), 

*  RHO( 20 ) . MTL( 20 ) . THETA( 20 ) . TH( 20 ) .QBAR( 5,6.20) 
C0MM0N/STF/ELXY(9,2),STIF(80.80),ELP(80).W0(80).W1(80),W2(80). 

*  A(5,5).B(3.3),D(3,3).AO,A1.A2.A3,A4,RH01.RH02.RH03 
COMMON/MSH/N0D(200.9).X(225),Y(225).DX(15).DY(15) 

DATA  NDF, NRMAX. NCMAX/5, 1000. 200/ 

C 

PI-3.14159265358 

C 

WRITE!*, '(427)' )  *  INPUT  THE  *.DAT  FILE  NAME  ’ 
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READ(*.*(A20)’)  OATFILE 

WRITE(*.'(A27)’)  •  TYPE  THE  OUTPUT  FILE  NAME  ‘ 

READ(*.*(A20)')  OUTPTF 
C 

lPRFIL-1 

IWRITE-2 

C 

OPEN  ( UNIT-I PRFI L. STATUS-’ OLD ’ . FORM-’ FORMATTED’ . FI LE-DATFI LE ) 
OPEN  (UNIT-IWRITE.STATUS-’NEM’ . FORM-’ FORMATTED ’, FI LE-OUTPTF) 

C 


PREPROCESSOR  UNIT 


Read  title  and  control  parameters  for  the  program 
READ(1.260)  TITLE 

READCl.*)  IEL.NPE.IMESH.NPRNT.ITEM.NTIME.NSTP.N02ER0 


General  Mesh  -  defined  by  element  connectivity  and  nodal  points 

IF  (IMESH.EQ.O)  THEN 
READd,*)  NEM.NNM 
00  10  I-l.NEM 

10  READd,*)  (NODd,J).J-l.NPE) 

READd,*)  (Xd),Yd),I-l,NNM) 

END  IF 

Rectangular  Generated  Mesh  -  generated  by  defining  number  of  element 
subdivisions  and  length  of  subdivisions  in  X  and  Y  directions 

20  IF  (IMESH.EQ.l)  THEN 
READd,*)  NX, NY 
NX1-IEL*NX 
NY1-IEL*NY 

READd,*)  (DXd).I-l,NXl) 

READd,*)  (DYd).I-l,NYl) 

CALL  MESHdEL,NX,NY,NPE,NNM,NEM) 

END  IF 

Read  the  properties  of  materials  used  in  the  plate 

30  READd,*)  NMTL 
DO  32  I-1,NMTL 

32  READd,*)  Eld),E2d),G12d),G13d),G23d).ANU12d),RH0d) 

Read  the  laminate  stacking  sequence  or  lay-up 

READd,*)  NLAY 
DO  34  I-1,NLAY 

34  READd,*)  MTLd),THETAd),THd) 

Read  pressure,  specified  forces,  and  specified  displacements 

READd.*)  PO 
READd,*)  NBDY 

READd.*)  dBDYd),I-l,NBDY) 

READd,*)  {VBDYd),I-l,NBDY) 
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READd,*)  NBSF 
IF  (NBSF.NE.O)  THEN 
READd,*)  (IBSFCn.I-l.NBSF) 
READd,*)  (VBSFd),I-l,NBSF) 
END  IF 


Transient  analysis  information  is  read 

3B  IF  (ITEM.EO.l)  THEN 
READd,*)  0T,ALFA 

Non-zero  initial  displacements  and  velocities 

IF  (NOZERO. EQ.l)  THEN 
NEQ-NNM*NDF 

READd,*)  (GFOd),I-l,NEQ) 

READd,*)  (GFld),I-l,NEQ) 

END  IF 

Time  Integration  parameters  (Newmark  scheme) 

36  BETA-0. 25* (0.5+ALFA)**2 
DT2-DT*0T 
A0-1.0/BETA/DT2 
A2-1.0/BETA/DT 
A1-ALFA*A2 
A3-0. 5/BETA- 1.0 
A4-ALFA/BETA-1.0 


Initialize  disp.,  vel .  and  accel .  vectors  if  not  specified 

IFCNOZERO.EQ.O)  THEN 
00  38  I-1,NEQ 
GF0d)-0.0 
GFld)-0.0 
GF2d)-0.0 
38  CONTINUE 
END  IF 
END  IF 


PROCESSOR  UNIT 


Compute  total  DOF's  'NEO',  and  element  DOF's  'NN' 

40  NE0-NNM*NDF 
NN-NPE*NDF 

Compute  the  plate  material  stiffness  and  inertial  properties 

CALL  MATPR0P{El.E2,G12,G13,G23,ANU12,-i.AY,MTL,THETA,RH0,TH.A,B,D, 
*  QBAR,H,RH01.RH02.RH03) 

Print  the  program  parameters  and  the  mesh  information 

WRITE(2,260)  TITLE 
WRITE(2,310)  IEL,NPE 


WRITE(2,320)  NEM.NNM.NOF 
C 

C  Print  the  material  properties  and  stacking  sequence 
C 

DO  50  I-1,NMTL 

50  MRITE(2.330)  I .El( I ) . E2( I) .612( I ) .G13( I ) ,G23( I ) . ANU12( 1 ) ,RH0( I ) 

WRITE(2,326) 

WRITE(2,327) 

DO  56  I-1,NLAY 

55  WRITE(2.328)  I .MTL( I ) .THETA( I )*180.0/PI .TH{ I ) 

WRITE(2.329)  H 
C 

C  Print  the  A,  B.  and  0  matrices  and  inertial  terms 
C 

WRITE(2,331) 

WRITE(2.332) 

WRITE(2.325)  ( A( 1 , J ) . J-1 ,3) 

WRITE(2.325)  (A(2. J ) . J-1 .3) 

WRITE{2.325)  (A(3.J) .J-1 .3) 

WRITE(2.334) 

WRITE(2,325)  (B(l . J ) , J-1 ,3) 

WRITE(2.325)  (B(2. J ) . J-1 .3) 

WRITE{2.325)  (B(3.J).J-1.3) 

MRITE(2.336) 

WRITE(2.325)  (0(1. J). J-1. 3) 

WRITE(2.325)  (0(2. J ) .J-1 .3) 

WRITE(2.325)  (0(3. J ) .J-1 .3) 

WRITE(2.338) 

WRITE(2.325)  A(4.4).A(4.5),A(5.5) 

WRITE(2.339) 

«R1TE(2.325)  RHOl . RH02 . RH03 
C 

C  Print  specified  displacements,  pressure  and  specified  forces 
C 

WRITE(2.345)  NBDY 

WRITE ( 2 . 280 )  ( I  BOY  a ) . 'B JY ( I ) . I-l . NBDY ) 

HRITE(2.342)  PO 
WR1TE(2.350)  NBSF 

WRITE(2.280)  ( IBSF( I ) . VBSF( I ) . I-l .NBSF) 

WRITE(2.360) 

C 

C  Print  element  connectivity  and  nodal  point  coordinates 
C 

DO  60  I-l.NEM 

60  WRITE(2.270)  I . (N0D( I . J) .J-1 .NPE) 

WRITE(2.370) 

WRITE(2.375)  ( I .X( I ) . Y( I ) .1-1 .NNM) 

C 

C  Compute  the  half-band  width  'NHBW'  'f  global  stiffness  matrix 
C 

NHBW-0 

DO  70  N-l.NEM 
00  70  I-l. NPE 
DO  70  J-1. NPE 

NW-( I ABS( N0D( N . I ) - N00( N . J ) )+l )*NOF 
IF  (NHBW.LT.NW)  NHBW-NN 
70  CONTINUE 

NRITE(2.400)NHBN 
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c 

T-0.0 

IF  (ITEM.EQ.l)  WRITE(2.460)  OT.ALFA.BETA.AO.Al .A2,A3.A4 
C 

C  .  00- Loop  on  number  of  time  steps  begins  here  . 

C 

DO  220  NT-1. NT I ME 

IF  (ITEM.EQ.l. AND. NT. GE.NSTP)  PO-0.0 


Initialize  the  global  stiffness  matrix  and  force  vector 

DO  80  I-l.NEQ 
GF(I)-0.0 
DO  80  J-l.NHBW 
80  GSTIF(I.J)-0.0 

Convert  global  information  to  the  element  level 

DO  130  N-l.NEM 
L-0 

DO  90  I-l.NPE 
NI-NOO(N.I) 

ELXY(I.1)-X(NI) 

ELXY(I.2)-Y(NI) 

LI-(NI-1)*NDF 
00  90  J-l.NDF 
LI-LI+1 
L-L+1 

W0(L)-GF0(LI) 

W1(L)-GF1(LI) 

W2(L)-GF2(LI) 

90  CONTINUE 

Compute  the  element  stiffness  and  mass  matrices 

CALL  STIFF! lEL.NPE.NN.PO. ITEM. NT. NOZERO) 

IF  (NPRNT.EQ.l)  THEN 
WRITE(2.380) 

00  100  I-l.NN 

100  WRITE  (2.300)  (STIF(I,J).J-1.NN) 

WRITE(2.410) 

WRITE(2.300)  (ELP(I). I-l.NN) 

WRITE(2.410) 

END  IF 


Assemble  element  stiffness  matrices  to  get  global  stiffness  matrix 

00  130  I-l.NPE 

NR-(N0D(N.I)-1)*NDF 
DO  130  II-l.NOF 
NR-NR+1 

L-(I-1)*NDF+II 

GF(NR)-GF(NR)+ELP(L) 

DO  130  J-l.NPE 
NCL-(N0D(N.J)-1)*NDF 
DO  130  JJ-l.NOF 
M-(J-1)*N0F+JJ 
NC-NCL+JJ-NR+1 


ooooo  r>r>o  ooo  r->or>o  oooo  or>or> 


IF  (NC.GT.O)  GSTIF(NR.NC)-6STIF(NR,NC)+STIF(L.M) 

130  CONTINUE 

The  global  system  equations  are  now  ready  for  implementing  the 
force  and  displacement  boundary  conditions 

IF  ((NBSF.GT.O). AND. (NOZERO. EQ.O. OR. ITEM. EQ.O))  THEN 
00  140  I-l.NBSF 
NB-IBSF(I) 

GF(NB)-GF(NB)+VBSF(I) 

140  CONTINUE 
END  IF 

145  CALL  BNDY(NRMAX.NCMAX,NEO.NHBN.GSTIF.GF.NBDY.IBDY,VBDY) 

Call  subroutine  SOLVE  to  solve  the  global  system  of  equations 
(the  solution  is  returned  in  GF(I)) 

CALL  SOLVE  (NRMAX.NCMAX.NEO.NHBW.GSTIF.GF.O) 

IF  (ITEM.EO.O)  GOTO  180 


Calculate  the  second  time  derivative  when  initial  conditions 
are  non-zero 

IF  ((NOZERO. EO.O). OR. (NT. GT.D)  GOTO  160 
00  150  I-l.NEO 
150  GF2(I)-GF(I) 

GOTO  210 

Calculate  new  velocities  and  accelerations 

160  T-T+DT 

00  170  I-INEO 

GF0( I )-A0*(GF( I ) -GF0( I ) ) -A2*GF1( I ) -A3*GF2( I ) 

GFl ( I )-GFl ( I )+DT*( 1 .0- ALFA )*GF2( I )+OT*ALFA*GFO(  I ) 
GF2(I)-GF0(I) 

GF0(I)-GF(I) 

170  CONTINUE 

Print  the  time  step  and  resulting  generalized  displacements 

WRITE(2,470)  T 
180  WRITE(2.480) 

WRITE(2.490)(((I+4)/NDF),GF(I).GF(I+l),GF(I+2),GF(I+3). 
*  GF(I+4).I-1.NE0.NDF) 

WRITE(2,410) 


POSTPROCESSOR  UNIT 


C  Compute  strains  and  stresses  (at  the  Gauss  points) 
C 

WRITE(2.440) 

WRITE(2.450) 

00  200  N-l.NEM 
L-0 

00  190  I-1,NPE 
NI-NOO(N,I) 


r>r>or>r)  ooo 


ELXY(I.1)-X(N1) 

ELXY(I.2)-Y(NI) 

LI-(N1-1)*NDF 
DO  190  J-l.NDF 
LI-LI+1 
L-L+1 

W0(L)-6F(LI) 

190  CONTINUE 

CALL  STRESS  (NPE.NOF.IEL.ELXY.NO.QBAR.NUY.TH.H) 
200  CONTINUE 

IF(ITEM.EQ.O)  GOTO  230 
210  WR1TE(2.390) 

220  CONTINUE 

.  End  of  DO-Loop  on  the  numtjr  of  time  steps  - 


230  STOP 


FORMATS 


260  FORMAT (20A4) 

270  F0RMAT(6X.I5.2X.9I5) 

280  F0RMAT((8X.3{2X.I4.2X.E12.5))) 

300  FORMAT(5(2X.E12.5)) 

310  FORMAT! /.5X/ ELEMENT  TYPE!  1-LINEAR. 2-QUADRATIC)  -’.I2.5X, 

*'  NODES  PER  ELEMENT-M2) 

320  FORMAT! lOX,' ACTUAL  NUMBER  OF  ELEMENTS  IN  THE  MESH-M3, 

*/.10X. 'NUMBER  OF  NODES  IN  THE  MESH  -'.13, 

*/.10X. 'DEGREES  OF  FREEDOM  -'.12./) 

325  FORMAT!7X.3!3X.E12.5)) 

326  FORMAT!/. 5X. 'LAMINATE  STACKING  SEQUENCE') 

327  FORMAT!/. 8X.' LAYER'. 3X.'MTL  #' .4X. 'THETA' .5X. 'THICKNESS' ) 

328  F0RMAT!8X.I3.6X.I2.2X.F10.5.2X.E12.5) 

329  FORMAT !/.8X. 'TOTAL  THICKNESS  -'.E12.5) 

330  FORMAT! 5X. 'MATERIAL  '.12.'  PROPERTIES: ' ./.lOX. 'MODULUS. El-' .E12. 5. 
*/.10X. 'MODULUS. E2-' .E12.5./.10X. 'SHEAR  MODULUS, G12-' .E12. 5. 

*/.10X. 'SHEAR  MODULUS. G13-'.E12. 5. /.lOX.' SHEAR  MODULUS. G23-' .E12.5. 
*/.10X. 'POISSONS  RATI0.NU12-' .E12.5. 

*/,10X. 'MATERIAL  DENSITY .RHO-' .E12. 5./) 

331  FORMAT!/. 5X. 'LAMINATE  PLATE  PROPERTIES') 

332  FORMAT !/.8X. 'A  MATRIX  TERMS') 

334  FORMAT!/. 8X.'B  MATRIX  TERMS') 

336  FORMAT! /.8X.'D  MATRIX  TERMS') 

338  FORMAT! /.8X. 'SHEAR  TERMS:  A44.  A45.  ASS') 

339  F0RMAT!/,8X.' INERTIAL  TERMS  RHOl.  RH02.  RH03') 

342  FORMAT! /.5X. 'UNIFORMLY  DISTRIBUTED  LOAD.  PO  -'.E12.5) 

345  FORMAT! /.5X. 'NUMBER  OF  SPECIFIED  DISPLACEMENTS-' . IS. 

♦/.SX.'SPECIFED  DISPLACEMENTS  AND  THEIR  VALUES  FOLLDW:') 

350  FORMAT! /.5X, 'NUMBER  OF  SPECIFIED  FORCES-' .14./. 5X. 

♦'SPECIFIED  FORCE  DEGREES  OF  FREEDOM  AND  THEIR  VALUES  FOLLOW:') 

360  FORMAT! /.5X. 'BOOLEAN  !CONNECTIVITY)  MATRIX-NOD!!. J)' ./) 

370  FORMAT!/. 5X. 'COORDINATES  OF  THE  GLOBAL  NODES:'./) 

375  F0RMAT!2!2X.I4.3X.E12.5,3X.E12.5)) 

380  FORMAT! /.5X. 'ELEMENT  STIFFNESS  AND  FORCE  MATRICES;'./) 

390  F0RMAT!120!' :').//) 


or>ooooooooooooor)ooor>oor>r>r>or>oo  ooooo 
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400  FORMAT (/.5X. ’HALF  BAND  WIDTH  OF  GLOBAL  STIFFNESS  MATRIX  -  MS,/) 
410  FORMAT!//) 

440  FORMATCSX,' LAMINATE  STRAINS  AND  STRESSES  AT  GAUSS  POINTS',/) 

450  FORMAT ( 1 X , ' ( X  - COORD , Y - COORD ) ' . 

*/,2X, *L0C' ,5X, 'EPSILONX* ,4X. 'EPSILONV .5X. 'GAMMAXY' .5X, 

* • KAPPAX • . 6X , • KAPPAY • , 6X , ■ KAPPAXY ' . / , 2X . ' LAY ' .4X , ' Z - COORD ' , 6X , 
**SIGMAX' ,6X,'SIGMAY' ,7X,'TAUXY',7X.’TAUY2’,7X.'TAUXZ'  ) 

460  FORMAT (/,5X,'DT-' .E10.4,5X. ' ALFA-' .E10.4,5X. 'BETA-' .E10.4,/ ,10X, 
‘•TEMPORAL  PARAMETERS  A0,A1 ,A2.A3.A4: ' ,5E12.4, / ) 

470  FORMAT!/, 5X, 'TIME-', ElO, 3./) 

480  FORMAT! /,5X, 'GENERALIZED  DISPLACEMENTS  !U.V,W.SX,SY)  PER  NODE',/) 
490  F0RMAT!2X,I4,2X,E12.5,2X.£12.5,2X.E12.5,2X,E12.5,2X,E12.5) 

END 

*********************************************************************** 

*  SUBROUTINES  * 

*********************************************************************** 


SUBROUTINE  MATPR0P!E1 ,E2,G12,G13. G23,ANU12.NLAY, MTL, THETA, RHO, 
T.A,B,0,0BAR,H.RH01.RH02,RH03) 


THIS  SUBROUTINE  CALCULATES  THE  0  AND  OBAR  MATRICES  FOR  EACH  LAYER 
AND  THE  LAMINATE  MATERIAL  PROPERTY  MATRICES  A!I.J).B!I ,J),D!I .J) 
AND  THE  INERTIAL  TERMS  RHOl .RH02.RH03 


A!I,J) . EXTENSIONAL  STIFFNESS  MATRIX  !I.J-1,2.3) 

A!K.L) . SHEAR  TERMS  OF  STIFFNESS  MATRIX  !K,L-4,6) 

AK . SHEAR  CORRECTION  COEFFICIENT 

AMM,ANN . SINE  AND  COSINE  OF  FIBER  ORIENTATION  'THETA' 

ANU12!I),ANU21!I),..P0ISS0N  RATIOS  OF  MATERIAL 

ATOL . ZERO  TOLERANCE  OF  STIFFNESS  TERMS  COMPARED  OTHERS 

B!I,J) . COUPLING  STIFFNESS  MATRIX  !I.J-1,2.3) 

D!I,J) . BENDING  STIFFNESS  MATRIX  !I,J-1.2,3) 

E1!I),E2!I).ELASTIC  MODULI  OF  MATERIAL  I 
612!I).G23!I),G13!I)..SHEAR  MODULI  OF  MATERIAL  I 

H . TOTAL  PLATE  THICKNESS 

MTL!L) . MATERIAL  NUMBER  OF  LAYER  L 

NLAY . NUMBER  OF  LAYERS  IN  THE  LAMINATE 

Q!I.J,L) _ STRESS-STRAIN  MARTIX  OF  UYER  L  ALIGNED  WITH 

PRINCIPAL  DIRECTIONS 

OBAR! I, J,L). TRANSFORMED  STRESS-STRAIN  MATRIX  OF  LAYER  L 

RHO! I) . DENSITY  OF  MATERIAL  I 

RHOl, RH02,RH03... INERTIAL  PARAMETERS  OF  THE  LAMINATE 

TH!L) . THICKNESS  OF  LAYER  L 

THETA! L).... FIBER  DIRECTION  ORIENTATION  OF  LAYER  L 
ZBAR!L) . MID-PLANE  POSITION  OF  LAYER  L 


IMPLICIT  REAL*8!A-H,0-Z) 

DIMENSION  Q!6,5,20).0BAR!5,5.20),E1!20).E2!20).G12!20),G13!20), 

*  G23!20) ,ANU12!20) ,ANU21!20) .THETA!20) ,T!20) , 

*  MTL!20),ZBAR!20).RHO!20).A(5.5).B!3.3),D!3,3) 

C 

C  Shear  correction  factor  'AK' 

C 


AK-5. 0/6.0 
H-0.0 


ooo  oor>  r>r>o 
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ATOL-l.OE-09 

PI-3.141592654 

C 

C  This  loop  calculates  'O'  &  'QBAR'  matrices  and  plate  thickness  'H' 

C 

DO  30  I-l.NLAY 
H-H+T(l) 

ANU21(MTL( I ) )-E2(MTL( I ) )*ANU12(MTL( I ) )/El(MTL( I ) ) 

DENOM-1 .0- ANU12(MTL( 1 ) )*ANU21 (MTU  I ) ) 

Q(1.1.I)-E1(MTL(I))/DEN0M 
Q(2.2.I)-E2(MTL(I))/0EN0M 
Q( 1 .2 . I )-ANU12(MTL( I ) )*E2(MTL( I ) )/0EN0M 
Q(3.3.I)-G12(MTL(n) 

Q(4.4.I)-G23(MTL(I)) 

Q(5.5.I)-G13(MTL(I)) 

Change  'THETA'  to  radians  and  calculate  cos  'AMM'  &  sin  'ANN' 

THETA( I )-THETA( I )*PI/180.0 
AMM-C0S(THETA(I)) 

ANN-SIN<THETA(n) 

Calculate  lamina  stiffness  matrices  'OBAR'  due  to  orientation  'THETA' 

QBAR(l.l.I)-Q(l.l.I)*(AMM*M.O)  +  2.0*(Q{1 .2. 1  )+2.0’‘Q(3.3. 1 )  )* 

*  {AMM*AMM)*(ANN*ANN)  +  Q(2.2.I )*(ANN**4.0) 
QBAR(1.2.I)-(0(1.1.I)+0(2.2.I)*4.0*(Q(3.3.I)))*(AMM*AMM)* 

*  (ANN* ANN )+  Q(1,2.I)*(AMM**4.0+ANN**4.0) 

QBARd. 3.1)— (0(2.2. I)*AMM*(ANN**3.0))  +  0(1 .1 . I )*(AMM**3.0)*ANN 

*  -  (0(l,2.I)+2.0*Q(3.3.I))*AWi*ANN*((AMM*AMM)-(ANN*ANN)) 
0BAR(2.1.1)-0BAR(1.2.I) 

OBAR(2.2.I)-0(1.1.I)*(ANN**4.0)  +  2.0*(0(1,2.I)+2.0*0(3.3.I))* 

*  (AMM*AMM)*(ANN*ANN)  +  0(2.2. I )*( AMM**4 .0) 

0BAR(2.3.I)— (0(2.2.I)*(AMM**3.0)*ANN)  +  0(1 .1 . 1 )*AMM*(ANN**3.0) 

*  +  (0(1.2.I)+2.0*0(3.3.I))*AMM*ANN*((AMM*AMM)-(ANN*ANN)) 
0BAR(3.1.I)-0BAR(1.3.I) 

0BAR(3.2.I)-0BAR(2.3.I) 

0BAR(3.3.I)-(0(1.1.I)+0(2.2.I)-2.0*0(1.2.I))*(AMM*AMM)*(ANN*ANN) 

*  +  0(3.3. I)*((AMM*AMM)-(ANN*ANN))**2.0 
QBAR(4.4.I)-Q(4.4.I)*(AMM*AMM)  +  0(5.5. I )*( ANN*ANN) 

0BAR(4 . 5 . I )-(0( 5 . 5 . I ) -0(4 .4 . I ) )*AMM*ANN 
0BAR(5.4.I)-QBAR(4.5.I) 

0BAR(5.5.I)-0(5.5.I)*(AMM*AMM)  +  0(4.4. I)*(ANN*ANN) 

30  CONTINUE 

Initialize  A.B.O  matrices  and  inertial  terms 

DO  40  1-1.3 
DO  40  J-1.3 
A(I.J)-0.0 
B(I.J)-0.0 
D(I.J)-0.0 
40  CONTINUE 
A(4.4)-0.0 
A(4.5)-0.0 
A(5.4)-0.0 
A(5.5)-0.0 
RHOl-0.0 


r>  o  o  o  r>  r> 


RH02-0.0 

RH03-0.0 


Calculate  lamina  mid-plane  position  'ZBAR'  from  laminate  mid-plane  (0) 

ZBAR(l)— H/2.0-H(l)/2.0 
IF(NLAY.GT.l)  THEN 
DO  50  I-2.NLAY 

ZBARC I  )-ZBAR(  I  - 1  )+T(  I  -1  )/2 .0-rT(  I  )/2 .0 
50  CONTINUE 
END  IF 

Calculate  A.B.D  in-plane  matrix  terms 

DO  80  K-1,NLAY 
00  60  1-1.3 
00  60  J-1.3 

A(I.J)-A(I.J)  -r  QBAR(I.J.K)*T{K) 

B(I.J)-B(I.J)  +  QBAR(I.J.K)*T(K)*ZBAR(K) 

D(I,J)-0(I.J)  QBAR(I.J.K)*{T(K)*2BAR(K)**2.0 

*  (T((0**3. 0/12.0)) 

60  CONTINUE 

Calculate  'A'  matrix  transverse  shear  terms 

DO  70  1-4.5 
DO  70  J-4.5 

A(I.J)-A(I.J)  +  AK*0BAR(I.J.K)*T(K) 

70  CONTINUE 

Calculate  inertial  terms  for  dynamic  analysis 

RHOl-RHOl  +  RHO(MTL(K))n(K) 

RH02-RH02  +  RH0(MTL(K) )*T(K)*ZBAR(K) 

RH03-RH03  +  RH0(MTL(K)  )*(T(K)*ZBAR(K)**2.0-KT(K)**3. 0/12.0) ) 

80  CONTINUE 

Set  negligibly  small  laminate  property  terms  to  zero 

T0L-A(1.1)*AT0L 
00  90  1-1.5 
00  90  J-1.5 

IF  (ABS(A(I,J)).LT.TOL)  A(I.J)-0.0 
IF  (I.LE.3.AND.J.LE.3)  THEN 
IF  (ABS(B(I.J)).LT.TOL)  B(I.J)-0.0 
IF  (ABS(D(I.J)).LT.TOL)  D(I.J)-0.0 
END  IF 
90  CONTINUE 

T0L-RH01*AT0L 

IF  (ABS(RH02).LT.T0L)  RH02-0.0 


RETURN 

END 


SUBROUTINE  ST I FF( I  EL . NPE . NN . PO . I TEM . NT . NOZERO ) 


Ill 

c 

c  . 

C  .  THIS  SUBROUTINE  IS  WRITTEN  FOR  COMPOSITE  PLATES.  THE  ELEMENT  IS  . 

C  .  BASED  ON  A  SHEAR-OEFORMABLE  THEORY.  HERE  THE  FOUR-.  EIGHT-  OR 
C  .  NINE-NODE  ISOPARAMETRIC  ELEMENT  WITH  FIVE  DEGREES  OF  FREEDOM 
C  .  (U.V.W.SX.SY)  PER  NODE  CAN  BE  USED  BY  SPECIFYING  THE  ElEMENT  TYPE  . 

C  . 

C  .  SUBROUTINE  VARIABLES 

C  . 

C  .  A(I,J) . EXTENSIONAL  STIFFNESS  MATRIX  (I. J-l, 2. 3) 

C  .  A(K,L) . SHEAR  TERMS  OF  STIFFNESS  MATRIX  {K.L-4.5) 

C  .  A0.A1.A2.A3.A4.. .PARAMETERS  IN  THE  TIME-APPROXIMATION  SCHEME 

C  .  B(I,J) . COUPLING  STIFFNESS  MATRIX  (I. J-l. 2, 3) 

C  .  CNST . INTEGRATION  CONSTANT  TRANSFORMED  TO  X.Y  COORDINATES  . 

C  .  D(I.J) . BENDING  STIFFNESS  MATRIX  (I. J-l. 2. 3) 

C  .  DET . DETERMINATE  OF  JACOBIAN  TRANSFORMATION  MATRIX 

C  .  ELP(I) . ELEMENT  FORCE  VECTOR 

C  .  ELXY( I. J)... ELEMENT  NODE  COORDINATES  OF  ELEMENT  NODE  I 
C  .  J-l  FOR  X-COORD.  J-2  FOR  Y-COORD 

C  .  GAUSS( I. J). .GAUSSIAN  POINT  COORDINATES  (LOCAL  VALUES) 

C  .  GOSF( J. I ).. .DERIVATIVE  OF  SHAPE  FUNCTION  SF( I) 

C  .  J-l  WITH  RESPECT  TO  X.  J-2  WITH  RESPECT  TO  Y 

C  .  H(I.J) . ELEMENT  MASS  MATRIX 

C  .  IBDY(I) . LOCATION  AND  DIRECTION  OF  SPECIFIED  GLOBAL 

C  .  DISPLACEMENTS 

C  .  IBSF(I) . LOCATION  AND  DIRECTION  OF  SPECIFIED  NONZERO  GLOBAL 

C  .  FORCES 

C  .  lEL . INDICATOR  FOR  THE  ELEMENT  TYPE: 

C  .  IEL-1.  4-NOOE  ELEMENT 

C  .  IEL-2,  8-  OR  9-NOOE  ELEMENT 

C  .  IMESH . INDICATOR  FOR  MESH  GENERATION 

C  .  (O-REAO  IN.  1-SOUARE  MESH  IS  GENERATED) 

C  .  ITEM . INDICATOR  FOR  TRANSIENT  ANALYSIS  (1-YES.  O-NO) 

C  .  LGP . ORDER  OF  REDUCED  INTEGRATION  ON  TRANSVERSE  SHEAR  TERM. 

C  .  NOF . NUMBER  OF  DEGREES  OF  FREEDOM  PER  NODE  (U.V  4.SX.SY)  . 

C  .  NGP . ORDER  OF  NORMAL  INTEGRATION  ON  IN-PLANE  TERMS 

C  .  NN . NUMBER  OF  DEGREES  OF  FREEDOM  PER  NODE 

C  .  (NODES  PER  ELEMENT  x  NODAL  DOF) 

C  .  NOZERO . INDICATOR  FOR  ZERO(NOZERO-O)  OR  NONZERO! NOZERO-1 ) 

C  .  INITIAL  CONDITIONS  FOR  TRANSIENT  ANALYSIS 

C  .  NPE . NUMBER  OF  NODES  PER  ELEMENT  (4.  8  OR  9) 

C  .  NT . CURRENT  TIME  STEP  NUMBER  IN  THE  TRANSIENT  ANALYSIS 

C  .  PO . INTENSITY  OF  APPLIED  TRANSVERSE  UNIFORM  PRESSURE 

C  .  SF(I) . VALUE  OF  INTERPOLATION  FUNCTION  OF  NODE  1 

C  .  STIF(I.J).. .ELEMENT  STIFFNESS  MATRIX 

C  .  SXX.SXY.SYY.SYX.. VALUES  OF  SHAPE  FUNCTION  DERIVATIVE  INTEGRALS 
C  .  SXO.SYO. SOX, SOY.. VALUES  OF  THE  PRODUCT  OF  SHAPE  FUNCTION  AND 
C  .  SHAPE  FUNCTION  DERIVATIVE  INTEGRALS 

C  .  SOO . VALUE  OF  SHAPE  FUNCTION  PRODUCT  INTEGRALS 

C  .  FOR  Snra  ABOVE  X-X  DERIVATIVE.  Y-Y  DERIVATIVE.  0-SHAPE  FUNCTION  . 

C  .  RH01,RH02,RH03.. LAMINATE  INERTIAL  PROPERTIES 

C  .  W0.N1.N2 _ ARRAYS  CORRESPONDING  TO  GFO.GF1.GF2  IN  AN  ELEMENT 

C  .  WT(I.J) . INTEGRATION  WEIGHT  VALUES 

C  .  XI. ETA . LOCAL  COORDINATE  VALUES  OF  GAUSS  POINTS 

C  . 

c 

IMPLICIT  REAL*8(A-H,0-Z) 

C0MM0N/STF/ELXY(9.2).STIF(80,80).ELP(80).W0(80).W1(80).W2(80). 
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*  A(5.5),B(3.3).D(3,3),A0.A1.A2.A3.A4.RH01,RH02.RH03 
C0MM0N/SHP/SF(9).GDSF(2,9) 

DIMENSION  GAUSS{4.4).WT(4.4).H(80.80) 

C 

C  Gaussian  Point  Coordinates 
C 

DATA  GAUSS/0. ODD,  0.000.  0.000.  0.000.  -.57735026918962600. 

*  .57735026918962600.  0.000.  0.000.  -.77459666924148300,  0.000, 

*  .77459666924148300.  0.000,  -.86113631159405300. 

*  -.33998104358485600.  .33998104358485600.  .86113631159405300/ 

C 

C  Integration  Weight  Values 
C 

DATA  WT/  2.000.  0.000.  0.000.  0.000,  1.000.  1.000.  0.000.  0.000. 

*  .55555555555555600,  .88888888888888900.  .55555555555555600. 

*  0.000.  .34785484513745400.  .65214515486254600. 

*  .65214515486254600,  .34785484513745400/ 

C 

C  Integration  order  of  in-plane  terms  ’NGP*.  and 
C  transverse  shear  terms  (reduced-integration)  ‘LGP' 

C 

NGP-IEL-H 

LGP-IEL 

NDF-NN/NPE 

C 

C  Initialize  the  element  matrices  'STIF'.  'H*  and  force  vector  'ELP' 

C 

00  10  1-1. NN 
ELP(I)-0.0 
00  10  J-1,NN 
H(I.J)-0.0 
STIF(I,J)-0.0 
10  CONTINUE 

Gauss  Quadrature  (Full  Integration)  on  in-plane  terms  begins  hare 

00  80  NI-l.NGP 
00  80  NJ-l.NGP 


Convert  to  local  Gauss  point  coordinates  and  evaluate  shape  function 

XI-GAUSS(NI.NGP) 

ETA-GAUSS(NJ,NGP) 

CALL  SHAPE(NPE,XI.ETA.ELXY,DET) 

CNST-DET*WT ( N I . NGP )*WT ( NO . NGP ) 

Distribution  of  constant  pressure  to  nodal  points 

00  30  I-1,NPE 
L-(I-l)*N0F-r3 

ELP(  L)-ELP(  L)-FCNST*SF(  I  )*P0 
30  CONTINUE 
C 

C  Compute  Stiffness  matrix  'STIF'  and  Mass  matrix  'H'  coefficients 
C 

II-l 

00  70  I-l.NPE 
JJ-1 


<-><-><->  o  o  o 


c 

c 

c 


DO  60  J-1,NPE 
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Integrals  of  the  shape  functions  and  derivatives 

SXX-GDSF(1 . 1 )*GDSF(1 .J)*CNST 
SYY-G0SF(2.I)*GDSF(2.J)*CNST 
SXY-GDSF(1.I)*GDSF{2.J)*CNST 
SYX-GDSF( 2 . I )*GDSF( 1 . J )*CNST 
S00-SF(I)*SF(J)*CNST 


Ful 1 -Integration  on  in-plane  stiffness  terms 

STIF(II.JJ)-STIF(II.JJ)+A{l,l)*SXX+A(1.3)*(SXY-rSYX)-HA(3.3)*SYY 

STIF(n.JJ+l)-STIF(Il,JJ+l)+A(1.2)«SXY+A(1.3)*SXX-HA(2,3)*SYY 

*  +A(3.3)*SYX 

STIF(II+l.JJ)-STIF(II+l.JJHA(l,2)*SYXrA(1.3)*SXX-rA(2,3)*SYY 

*  +A(3,3)*SXY 

ST  I F  ( 1 1 .  J  J-r3  )-ST  I F  ( 1 1 .  J  J+3  )+B  ( 1 , 1 )  *SXX-t-B  ( 1 . 3 )  *  ( SXY+SY  X ) 

*  +B(3.3)*SYY 

STI  F(  I  I-r3 .  JJ  )-STI  F(  1 1+3 .  JJ  )+B(  1 . 1  )*SXX+B(  1 . 3  )*(  SYX+SXY ) 

*  +B(3,3)*SYY 

STIF(II,JJ+4)-STIF(II.JJ+4)+B(l,2)*SXY+B(1.3)*SXX+B(2.3)*SYY 

*  +B(3.3)*SYX 
STIF(II+4,JJ)-STIF(II+4.JJ)+B(1.2)*SYX+B(1,3)*SXX+B(2.3)*SYY 

+  +B(3.3)*SXY 

STIF(II+1.JJ+1 )-STIF(II+l.JJ+l)+A(2.2)*SYY+A(2.3)*(SXY+SYX) 

+  +A(3.3)*SXX 

STIF(II+l.JJ+3)-STIF(II+l.JJ+3)+B(l,2)*SYX+B(2.3)*SYY 

*  +B(1.3)*SXX+B(3.3)*SXY 
STIF(II+3.JJ+1 )-STIF(II+3.JJ+l )+B(l .2)*SXY+B(2.3)*SYY 

*  +B(1.3)*SXX+B(3.3)*SYX 
STIF(II+l.JJ+4)-STIF(II+l,JJ+4)+B(2.2)*SYY+B{2.3)*(SXY+SYX) 

*  +B(3.3)*SXX 
STIF(II+4,JJ+1)-STIF(II+4,JJ+1)+B(2.2)*SYY+B(2.3)*(SYX+SXY) 

*  +B(3.3)*SXX 
STIF(n+3.JJ+3)-STIF(II+3.JJ+3)+D(l.l)*SXX+D(l,3)*(SXY+SYX) 

*  +0(3.3)*SYY 

STI F( 1 1+3 . JJ+4 )-STI F( 1 1+3 . JJ+4 )+0( 1 .2 )*SXY+D( 1 , 3 )*SXX 

*  +0(2.3)*SYY+0(3,3)*SYX 
STIF(II+4,JJ+3)-STIF(II+4,JJ+3)+0(1.2)*SYX+D(1.3)*SXX 

*  +D{2.3)*SYY+0(3,3)*SXY 
STIF(II+4,JJ+4)-STIF(II+4.JJ+4)+0(2,3)*(SXY+SYX)+D(3.3)*SXX 

*  +D(2.2)*SYY 


Mass  matrix  terms  'H'  for  transient  analysis 

IF  (ITEM. £0.1)  THEN 
H( 1 1 , JJ )-H( 1 1 , JJ )+RH01*S00 
H{ 1 1 . JJ+3 )-H( 1 1 , JJ+3 )+RH02*S00 
H( 1 1+3 . JJ )-H( 1 1+3 .JJ )+RH02*S00 
H( I I+l , JJ+I )-H( I I+l . JJ+1 )+RH01*SOO 
H{ I I+l , JJ+4 )-H( I I+l . JJ+4)+RHO2*S00 
H( 1 1+4 , JJ+1 )-H( 1 1+4 .JJ+1 )+RH02*S00 
H( 1 1+2 . J J+2 )-H( 1 1+2 . JJ+2 )+RH01*S00 
H( 1 1+3 .JJ+3 )-H( 1 1+3 .JJ+3 )+RH03*S00 
H( 1 1+4 .JJ+4 )-H( 1 1+4 .JJ+4 )+RH03+S00 
END  IF 


ooo  o  or>o  r>r>r>  r>  r>oo 
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JJ-NDF*J+1 
60  CONTINUE 

II-NDF*I+1 
70  CONTINUE 
80  CONTINUE 

Gauss  Quadrature  (Reduced  Integration)  on  transverse  shear  terms 

DO  no  NI-1.L6P 
DO  110  NJ-l.LGP 
XI-GAUSS(NI.L6P) 

ETA-GAUSS(NJ.L6P) 

CALL  SHAPEINPE.Xl.ETA.ELXY.OET) 

CNST-DET*WT( NI . LGP )*WT( NJ . L6P ) 

II-l 

DO  100  I-1,NPE 
JJ-1 

DO  90  J-l.NPE 

Integrals  of  the  shape  functions  and  derivatives 

SXX-GDSF( 1 . I )*G0SF( 1 . J )*CNST 

SYY-GDSF(2 . I )*G0SF( 2 , J )*CNST 

SXY-GDSF(1,I)*GDSF(2.J)*CNST 

SYX-GDSF( 2 . I )*GDSF( 1 . J )*CNST 

SX0-GDSF(1.I)*SF(J)*CNST 

S0X-SF(I)*GDSF(1.J)*CNST 

SYO-GOSF(2.I)*SF{J)*CNST 

S0Y-SF(I)*G0SF(2,J)*CNST 

S00-SF(I)*SF(J)*CNST 

Reduced- Integration  on  in-plane  stiffness  terms 

STIF(II+2.JJ+2)-STIF(II-r2.J0-r2)-t-A{5.5)*SXX+A(4,5)*(SXY+SYX) 
*  +A(4.4)*SYY 

STIF(II+2,JJ+3)-STIF(II+2,JJ-^3)+A(6.5)*SX0+A(4,5)*SY0 
STIF(II+3.JJ+2)-STIF(II+3.JJ+2)+A{5.6)*S0X+A(4.5)*S0Y 
STI F( 1 1+2 , JJ+4 )-STI F( 1 1+2 , JJ+4 )+A(4 , 5 )*SX0+A(4 .4 )*SY0 
STIF( 1 1+4 . JJ+2 )-STIF( 1 1+4 . JJ+2)+A(4 ,5)*S0X+A(4 ,4 )*S0Y 
ST I F( 1 1+3 . J J+3 )-ST I F( 1 1+3 . J J+3 )+A( 5 . 5 )*S00 
ST I F( 1 1+3 . JJ+4 )-ST I F( 1 1+3 .JJ+4 )+A(4 , 5 )*S00 
ST I F( I 1+4 . J J+3 )-ST I F( I 1+4 , J J+3 )+A( 4 . 5 )*S00 
STI F( 1 1+4 , JJ+4 )-STI F( 1 1+4 , JJ+4 )+A(4 .4 )*S00 

JJ-NDF*J+1 
90  CONTINUE 

II-NDF+I+1 
100  CONTINUE 
110  CONTINUE 


Element  calculations  for  transient  analysis  brgin  here 

IF  (ITEM.EQ.O)  RETURN 
IF  (NOZERO. EQ.l. AND. NT. EQ.l)  THEN 
DO  120  I-l.NN 
ELP(I)-0.0 
DO  120  J-l.NN 


ELP( I )-ELP( I )-STIF{ I . J)*WO(J) 

STIF(I.J)-H(I.J) 

120  CONTINUE 
RETURN 
END  IF 

130  00  140  I-l.NN 
00  140  J-l.NN 

ELP( I )-ELP( I )+H( I , J )*( AO*WO{ J )+A2*Wl ( J )+A3*W2(  J ) ) 
STIF(I.J)-STIF(I.J)+AO*H(I.J) 

140  CONTINUE 
RETURN 
END 


k ************************ 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  STRESS  (NPE.NDF.IEL.ELXY.W.QBAR.NLAY.TH.H) 


THIS  ROUTINE  EVALUATES  THE  STRESSES  AND  STRAINS  AT  THE  GAUSS 
POINTS  USING  THE  REDUCED  INTEGRATION. 

AKAPX.AKAPY.AKAPXY.. CURVATURES  AT  CURRENT  GAUSS  POINT 
ELXY(I.J)... ELEMENT  NODE  COORDINATES  OF  ELEMENT  NODE  I 
J-1  FOR  X- COORD.  J-2  FOR  Y- COORD 

EPN(I) . VECTOR  OF  CURRENT  2-POSITION  STRAINS 

l-EPNxx.2-EPNyy,3-GAMxy.4-GAMyz.6-GAMxz 
EPNXO.EPNYO.GAMXYO.. MID-PLANE  STRAINS  AT  CURRENT  GAUSS  POINT 
EPNXl.EPNYl.GAMXYl.. LAMINATE  BOTTOM  STAINS  AT  CURRENT  GAUSS  POINT 
EPNX2.EPNY2.GAM^Y2.. LAMINATE  TOP  STAINS  AT  CURRENT  GAUSS  POINT 
GAMYZ.GAMXZ,. TRANSVERSE  SHEAR  STRAINS  AT  CURRENT  GAUSS  POINT 


GAUSS( I. J).. GAUSSIAN  POINT  COORDINATES  (LOCAL  VALUES) 

GDSF(J. I)... DERIVATIVE  OF  SHAPE  FUNCTION  SF(I) 

J-1  WITH  RESPECT  TO  X,  J-2  WITH  RESPECT  TO  Y 

H . TOTAL  PLATE  THICKNESS 

lEL . INDICATOR  FOR  THE  ELEMENT  TYPE: 

I  EL-1.  4 -NODE  ELEMENT 
IEL-2,  8-  OR  9-NOOE  ELEMENT 

L . POINTER  TO  FIRST  OGF  OF  NODE  WITH  'NDF'  OOF  PER  NODE 

NOF . NUMBER  OF  DEGREES  OF  FREEDOM  PER  NODE  (U.V.W.SX.SY) 

NGP . ORDER  OF  REDUCED  INTEGRATION  FOR  STRAINS 

NLAY . NUMBER  OF  PLATE  LAYERS 

NPE . NUMBER  OF  NODES  PER  ELEMENT  (4.  8  OR  9) 

OBAR( I. J.L). TRANSFORMED  STRESS-STRAIN  MATRIX  OF  LAYER  L 

SIGMA(I) _ VECTOR  OF  CURRENT  Z-POSITION  STRESSES 

1-Sxx.2-Syy ,3-Txy.4-Tyz.6-Txz 

SF(I) . VALUE  OF  INTERPOUTION  FUNCTION  OF  NODE  I 

TH(L) . THICKNESS  OF  LAYER  L 

W(I) . VALUES  OF  ELEMENT  GENERALIZED  DISPLACEMENTS 

XI. ETA . LOCAL  COORDINATE  VALUES  OF  GAUSS  POINTS 

X.Y . GLOBAL  COORDINATES  OF  CURRENT  GAUSS  POINT 

Z(L) . Z-COORDINATE  OF  LAYER  INTERFACES  (LAYER  L/2) 

L-ODD-BOTTOM  OF  LAYER.  L-EVEN-TOP  OF  LAYER 


IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON/SHP/SF( 9 ) . GOSF( 2 , 9 ) 

DIMENSION  GAUSS(4.4).ELXY(9.2).W(80).QBAR(5.6.20).TH(20).Z(40), 


ooo  ooo  ooor> 


EPN(5).SIGMA{5) 


* 


C  Gaussian  point  coordinates 
C 

DATA  GAUSS/0.000.  O.ODO.  0.000.  0.000.  -.57735026918962600. 

*  .57735026918962600.  0.000.  0.000.  -.77459666924148300.  0.000. 

*  .77459666924148300.  O.ODO.  -.86113631159405300, 

*  -.33998104358485600.  .33998104358485600.  .86113631159405300/ 

C 

C  Order  of  integration  (reduced)  'NGP* 

NGP-IEL 

00  60  NI-1,NGP 
00  60  NJ-l.NGP 
XI-GAUSS(NI.NGP) 

ETA-GAUSS(NJ.N6P) 

CALL  SHAPE  (NPE.XI.ETA.ELXY.OET) 

EPNXO-0.0 
EPNYO-0.0 
GAMXYO-0.0 
GAMY 2-0.0 
GAMX2-0.0 
AKAPX-0.0 
AKAPY-0.0 
AKAPXY-0.0 
X-0.0 
Y-0.0 

DO  20  I-l.NPE 
L-(I-1)*N0F+1 
X-X+SF(I)*ELXY(I.l) 

Y-Y+SF(n*ELXY(I.2) 

Compute  the  midplane  strains  and  curvatures  and  average 
tranverse  shear  strain 

10  EPNX0-EPNX0+G0SF(1.I)*H(L) 

EPNYO-EPNYO+GOSF(2. I )*W( L+1 ) 

GAMXY0-GAMXY0+GDSF(2 , I )*W( L)+GOSF( 1 , 1 )*W( L+1 ) 
GAMYZ-GAMYZ+SF( I )*W( L+4 )+6DSF(2 . I )+W( L+2 ) 

GAMXZ-GAMXZ+SF ( I )*W( L+3)+6DSF( 1 . I )*M( L+2 ) 
AKAPX-AKAPX+GDSF( 1 . I )*W( L+3 ) 

AKAPY-AKAPY+GDSF(2.I )*W(L+4) 

AKAPXY-AKAPXY+GOSF( 2 , I )*H( L+3 )+GDSF{ 1 , I )*W( L+4 ) 

20  CONTINUE 

Compute  strains  at  the  bottom  (1)  and  top  (2)  of  the  laminate 

EPNX1-EPNX0-(H/2.0)*AKAPX 
EPNY1-EPNY0-(H/2.0)*AKAPY 
GAMXY1-6AMXY0-(H/2.0)*AKAPXY 
EPNX2-EPNX0+( H/2 .0 )*AKAPX 
EPNY2-EPNY0+( H/2 . 0 ) * AKAPY 
GAMXY2-GAMXY0+(H/2 . 0 )*AKAPXY 

Print  midplane  strains  and  curvatures,  and  surface  strains 
WRITE(2,100)  X.Y 

WRITE(2.110)  EPNXO.EPNYO.GAMXYO.AKAPX.AKAPY.AKAPXY 


ooooooooooo  oo  ooo  ooo 


WRITE(2.120)  EPNXl.EPNYl.GAMXYl 
WRITE(2.130)  EPNX2.EPNY2.GAMXY2 
C 

C  Compute  z -coordinate  at  the  lamina  interfaces 
C 

Z(0)— H/2.0 
00  50  LL-1,NLAY 
2(2*LL-l)-2(2*LL-2) 

Z(2*LL)-Z(2*LL-1)+TH(LL) 

C 

C  Compute  layer  strains  at  layer  interfaces 
C 

00  40  KK-1.2 

EPN(  1  )-EPNX0-rZ{2*(  LL-1  )+KK)*AKAPX 
EPN(2)-EPNY0+Z(2*(LL-l)+KK)*AKAPy 
EPN(  3  )-GAMXY0-t-Z(  2*(  LL- 1  )-H(K)*AKAPXY 
EPN(4)-GAHYZ 
EPN(5)-GAMX2 

This  loop  computes  the  layer  stresses  at  layer  interfaces 

00  30  11-1,5 
SlGMA(II)-0.0 
00  30  JJ-1.5 

SIGMA( 1 1 )-SlGMA( 1 1 )+QBAR( 1 1 . JJ , LL)*EPN( JJ ) 

30  CONTINUE 

Print  the  stresses  at  the  lamina  interfaces 

WRITE(2.140)  LL,Z(2*(LL-1)+KK).(SIGMA(MM),MM-1.5) 
40  CONTINUE 

50  CONTINUE 
60  CONTINUE 
RETURN 


100  F0RMAT(/,1X.'C  ,E12.«.  ’  .•.E12.4.')’) 
110  F0RMAT(2X,'MIDMX.riX,.SE12.4)) 

120  F0RMAT(2X.  'BOT'  .1X.(1/..3E12.4)) 

130  F0RMAT(2X,’T0PMX.(1X,3E12.4)) 

140  F0RMAT(2X,I2.1X.E12.4.1X.5E12.4) 

ENO 


SUBROUTINE  BNOY ( NRMAX , NCMAX , NEQ , NHBM . S . SL. NBOY , I BOY . VBOY ) 


SUBROUTINE  USEO  TO  IMPOSE  BOUNOARY  CONDITIONS  ON  BANOED  EQUATIONS 


IBOY(I) . LOCATION  ANO  OIRECTION  OF  SPECIFIED  GLOBAL 

DISPLACEMENTS 

NBDY . TOTAL  NUMBER  OF  SPECIFIED  GLOBAL  DISPLACEMENTS 

NCMAX . VALUE  OF  THE  COLUMN-DIMENSION  OF  S 

C  .  NEQ . TOTAL  NUMBER  OF  DEGREES  OF  FREEDOM  (NOOESxNODAL  DOF) 

C  .  NHBW . HALF-BAND  WIDTH  OF  GLOBAL  STIFFNESS  MATRIX 

C  .  NRMAX . VALUE  OF  THE  ROW-DIMENSION  OF  S 

C  .  S(M,N) . GLOBAL  STIFFNESS  MATRIX  (IN  BANDED  FORM) 


ooooooooooooooooooo  o  ooooo 


.  SUM) . GLOBAL  FORCE  VECTOR 

.  SVAL . VALUE  OF  CURRENT  SPECIFIED  DISPLACEMENT 

.  VBDY(I) . VALUES  OF  THE  DISPLACEMENTS  CORRESPONDING  TO  IBDYCI)  . 


IMPLICIT  REAL*8(A-H.0-Z) 

DIMENSION  S( NRMAX . NCMAX ) . SL( NRMAX ) . I BOY ( NBDY ) . VBOY ( NBDY ) 

DO  300  NB-l.NBOY 
lE-IBDY(NB) 

SVAL-VBDY(NB) 

IT-NHBW-1 
I-IE-NHBW 
DO  100  II-l.IT 
I-I+I 

IF  (I.GE.l)  THEN 
J-IE-I+1 

SL(I)-SL(I)-S(I.J)*SVAL 
S(I.J)-0.0 
END  IF 
100  CONTINUE 
S(I£.I)-1.0 
SL(IE)-SVAL 
I-IE 

00  200  II-2.NHBN 
I-I+l 

IF  (I.LE.NEQ)  THEN 
SL(I)-SL(I)-S(IE.II)*SVAL 
S(IE.Il)-0.0 
END  IF 
200  CONTINUE 
300  CONTINUE 
RETURN 
END 


SUBROUTINE  SOLVE ( NRM . NCM . NEONS . NBW. BAND. RHS , I  RES  > 


SOLVING  A  BANDED  SYMMETRIC  SYSTEM  OF  EQUATIONS 
IN  RESOLVING,  IRES  .GT.  0.  LHS  ELIMINATION  IS  SKIPPED 

BANO(M,N).,. GLOBAL  STIFFNESS  MATRIX  (IN  BANDED  FORM) 


IRES . IF  IRES  .GT.  0  THEN  FORWARD  ELIMINATION  IS  SKIPPED 

NBW . HALF-BAND  WIDTH  OF  GLOBAL  STIFFNESS  MATRIX 

NCM . VALUE  OF  THE  COLUMN-DIMENSION  OF  S 

NEONS . NUMBER  OF  EQUATIONS  (TOTAL  DEGREES  OF  FREEDOM) 

NRM . VALUE  OF  THE  ROW-OIMENSION  OF  S 

RHS(M) . GLOBAL  FORCE  VECTOR 

SVAL . VALUE  OF  CURRENT  SPECIFIED  DISPLACEMENT 

VBDY(I) . VALUES  OF  THE  DISPLACEMENTS  CORRESPONDING  TO  IBOY(I) 


IMPLICIT  REAL*8(A-H.0-Z) 
DIMENSION  BAND(NRM.NCM),RHS(NRM) 
C 

MEQNS-NEQNS-1 


or^o  r>r>o  ooo 
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C 

C  Forward  Elimination 
C 

IF  (IRES.EQ.O)  THEN 
00  500  NPIV-l.MEQNS 
NPIVOT-NPIV+1 
LSTSUB-NPIV+NBM-l 
IF  (LSTSUB.GT. NEONS)  LSTSUB-NEQNS 
00  400  NROM-NPIVOT.LSTSUB 


Invert  rows  and  columns  for  row  factor 
NCOL-NROH-NPIV+1 

FACTOR-B AN0( NP I V . NCOL ) / B AN0( NP I V . 1 ) 

00  200  NCOL-NROW.LSTSUB 
ICOL-NCOL-NROW+1 
JCOL-NCOL-NPIV+1 

BANO(NROW.ICOL)-BANO(NRON.1COL)-FACTOR*BANO(NPIV.JCOL) 
200  CONTINUE 

RHS ( NROW )-RHS ( NROW ) - FACTOR*RHS( NP I V ) 

400  CONTINUE 
500  CONTINUE 
ENO  IF 


Skip  forward  elimination  of  matrix  for  IRES  .GT.  0 

IF  (IRES.GT.O)  THEN 
90  00  100  NPIV-l.HEONS 

NPIVOT-NPIV+1 
LSTSUB-NPIV+NBW-1 
IFdSTSUB.GT. NEONS)  LSTSUB-NEQNS 
00  110  NROW-NPIVOT.LSTSUB 
NCOL-NROW-NPIV+1 

FACT0R-BAN0( NP I V . NCOL) /BAN0( NPI V . 1 ) 
RHS(NROH)-RHS(NROW)-FACTOR*RHS(NPIV) 

110  CONTINUE 
100  CONTINUE 
ENO  IF 


Back  Substitution 

00  800  IJK-2.NE0NS 
NPIV-NEQNS-IJK+2 
RHS(NPIV)-RHS(NPIV)/BANO{NPIV.l) 
LSTSUB-NPIV-NBW+1 
IF  (LSTSUB.LT.l)  LSTSUB-1 
NPIVOT-NPIV-1 
00  700  JKI-LSTSUB.NPIVOT 
NROW-NPIVOT-JKI+LSTSUB 
NCOL-NPIV-NROH+1 
FACT0R-BAN0( NROW .NCOL) 
RHS(NROW)-RHS(NROH)-FACTOR*RHS{NPIV) 
700  CONTINUE 
800  CONTINUE 

RHS(1)-RHS(1)/BAN0(1,1) 

RETURN 

END 


ooo  o  r>r)r>  or>o  r^oooooooooooo 


Q************************************************************************ 

c 

SUBROUTINE  SHAPECNPE.XI .ETA.ELXY.DET) 

C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 


THIS  SUBROUTINE  EVALUATES  THE  INTERPOLATION  FUNCTIONS  (SFd))  AND 
ITS  DERIVATIVES  WITH  RESPECT  TO  NATURAL  COORDINATES  (DSFCI.J)). 
AND  THE  DERIVATIVES  OF  SFd)  WITH  RESPECT  TO  GLOBAL  COORDINATES 
FOR  4.  8.  AND  9-NODED  RECTANGULAR  ISOPARAHETIRC  ELEMENTS 

DET . DETERMINATE  OF  JACOBIAN  TRANSFORMATION  MATRIX 

DSFCI.J) _ LOCAL  DERIVATIVE  OF  SF(J)  WITH  RESPECT  TO  XI  IF  I-l 

AND  WITH  RESPECT  TO  ETA  IF  1-2. 

ELXYCI.J)... ELEMENT  NODE  COORDINATES  OF  ELEMENT  NODE  I 
J-1  FOR  X-COORD,  J-2  FOR  Y-COORD 
GDSFCI.J)... GLOBAL  DERIVATIVE  OF  SF(J)  WITH  RESPECT  TO  X  IF  I-l 
AND  WITH  RESPECT  TO  Y  IF  1-2. 


GJCI.J) . JACOBIAN  MATRIX 

GIINVCI.J).. INVERSE  OF  THE  JACOBIAN  MATRIX 

NPCI) . ARRAY  OF  ELEMENT  NODES  (USED  FOR  DEFINING  SF  AND  DSF) 

NPE . NUMBER  OF  NODES  PER  ELEMENT  (4.  8  OR  9) 

SFd) . INTERPOLATION  FUNCTION  FOR  NODE  I  OF  THE  ELEMENT 

XI. ETA . LOCAL  COORDINATE  VALUES  OF  GAUSS  POINTS 

XNOOEd.J)..LOCAL  COORDINATES  OF  NODE  I  OF  THE  ELEMENT 
J-1  FOR  XI -COORD,  J-2  FOR  ETA- COORD 


IMPLICIT  REAL*8  (A-H.O-Z) 

C0MM0N/SHP/SF(9).GDSF(2.9) 

DIMENSION  ELXY(9.2).XN0DE(9.2).NP(9).DSF(2.9).GJ(2.2).GJINV(2.2) 

Local  nodal  point  coordinates  'XNOOE'  and  node  numbers  ’NP’ 

DATA  XNODE/-1. ODO. 2*1. ODO. -1.000.0. ODD. 1. ODD. 0. OOO. -1. ODD. O.ODO, 
*  2*-l. ODO. 2*1. ODO. -1.000. 0.000. 1.000. 2*0. ODO/ 

DATA  NP/1.2,3.4,5,7.6,8.9/ 

Multiplication  function  for  real  variables 

FNC(A,B)-A*B 

IF  (NPE-8)  60.10.80 

Quadratic  interpolation  functions  (for  the  EIGHT-NODE  element) 

10  DO  40  I-l. NPE 
NI-NPCI) 

XP-XNODE(NI.l) 

YP-XN0DE(NI.2) 

XI0-1.0+XI*XP 

ETA0-1.0+ETA*YP 

XI1-1.0-XI*XI 

ETA1-1.0-ETA*ETA 

IF  (I.GT.4)  GOTO  20 

SF(NI)-0.25*FNC(XI0,ETA0)*(XI*XP+ETA*YP-1.0) 

DSF(1.NI)-0.25*FNC(ETAO.XP)*(2.0*XI*XP+ETA*YP) 

DSF(2.NI)-0.25*FNC(XI0.YP)*(2.0*ETA*YP+XI*XP) 


O  O  O  O  O  <-) 


I2I 

GOTO  40 

20  IF  (I.GT.6)  GOTO  30 

SF{NI)-0.5*FNC{Xn.ETA0) 

DSF(l.NI)— FNC(Xl.ETAO) 

DSF(2.NI)-0.5*FNC(YP.XI1) 

GOTO  40 

30  SF{NI)-0.5*FNC(ETA1.XI0) 

0SF(1.NI)-0.5*FNC(XP.ETA1) 

0SF(2.NI)— FNC(ETA.XIO) 

40  CONTINUE 
GOTO  130 

Linear  interpolation  functions  (for  the  FOUR-NODE  element) 

60  DO  70  I-l.NPE 
XP-XNODEd.l) 

YP-XN00E(I,2) 

XI0-1.0+XI*XP 

ETA0-1.0+ETA*YP 

SF(I)-O.25*FNC(XI0.ETA0) 

0SF(1.I)-0.25*FNC(XP.ETA0) 

0SF(2.I)-0.25*FNC(YP,XI0) 

70  CONTINUE 
GOTO  130 

Quadratic  interpolation  functions  (for  the  NINE-NODE  element) 

80  DO  120  I-l.NPE 
NI-NP(I) 

XP-XNOOE(NI.l) 

YP-XN00E(NI.2) 

XI0-1.0+XI*XP 
ETA0-1.0+ETA*YP 
XI1-1.0-XI*XI 
ETA1-1.0-ETA*ETA 
XI2-XP*XI 
ETA2-YP*ETA 
IF  (I.GT.4)  GOTO  90 
SF(NI)-0.25*FNC(XI0.ETA0)*XI2*ETA2 
0SF( 1 , N I )-0 . 25*XP*FNC( ETA2 . ETAO )*( 1 . 0+2 . 0*X I 2 ) 
DSF(2.NI)-0.25*YP*FNC(XI2,XIO)*(1.0+2.0*ETA2) 

GOTO  120 

90  IF  (I.GT.6)  GOTO  100 

SF(NI)-0.5*FNC(XI1.ETA0)*ETA2 
OSF(l.Nl)— XI*FNC(ETA2.ETA0) 

DSF(2,NI)-0.5+FNC(XI1,YP)*(1.0+2.0*ETA2) 

GOTO  120 

100  IF  (I.GT.8)  GOTO  110 

SF(NI)-0.5*FNC(ETA1.XI0)*XI2 
0SF(2.NI)— ETA*FNC(XI2.XI0) 

DSFd  ,NI  )-0.5*FNC( ETAl  ,XP)*(1 .0+2.0+XI2) 

GOTO  120 

110  SF(NI)-FNC(XI1,ETA1) 

DSFd.ND— 2.0+XI+ETAl 
DSF(2.NI)— 2.0*ETA*XI1 

120  CONTINUE 


o  o  o 


c 

c 

c 


Transform  derivatives  from  local  (XI. ETA)  to  global 


130 


140 


150 


DO  140  1-1.2 
DO  140  J-1.2 
GJ{I.J)-0.0 
DO  140  K-l.NPE 

GJ(I.J)-6J(I.J)+DSF(I.K)*ELXY(K.J) 

CONTINUE 

0ET-GJ(1.1)*GJ(2.2)-GJ(1.2)*GJ(2.1) 
GJINV(1.1)-GJ(2.2)/0ET 
GJINV(2.2)-GJ(1.1)/0ET 
GJINV(1.2)— GJ(1.2)/0ET 
GJINV(2.1)— GJ(2.1)/DET 
DO  150  1-1.2 
DO  150  J-l.NPE 
G0SF(I.J)-0.0 
DO  150  K-1.2 

GOSF(I.J)-GOSF(I.J)+GJINV(I.K)*DSF(K.J) 

CONTINUE 


(X.Y)  derivatives 


C 

RETURN 

END 

C 


C 

SUBROUTINE  MESH(IEL.NX.NY.NPE.NNM.NEM) 

C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON/MSH/NOO(200.9).X(225).Y(225).DX(15).DY(15) 


THIS  SUBROUTINE  GENERATES  ARRAY  NOOCI.J)  COORDINATES  X(I).Y(I) 

AND  MESH  INFORMATION  (NNM.NEM.NPE)  FOR  RECTANGULAR  DOMAINS.  THE 
DOMAIN  IS  DIVIDED  INTO  LINEAR  OR  QUADRATIC  QUADRILATERAL  ELEMENTS 
(NX  BY  NY  NONUNIFORM  MESH  IN  GENERAL). 


DX( I ).0Y( I). DISTANCE  BETWEEN  NODES  IN  X.Y  DIRECTIONS  FOR  MESH 
GENERATION 

lEL . ELEMENT  TYPE  (IEL-1:  4  NODES.  IEL-2:  8  OR  9  NODES) 

NEM . TOTAL  NUMBER  OF  ELEMENTS 

NNM . TOTAL  NUMBER  OF  NODES 

NOD( I. J)... CONNECTIVITY  MATRIX 

NPE . NUMBER  OF  NODES  PER  ELEMENT 

NX. NY . NUMBER  OF  ELEMENTS  ALONG  X.Y-DIRECTIONS 

NXX.NYY _ NUMBER  OF  SUBDIVISIONS  BETWEEN  NODES  IN  X.Y-DIRECTIONS 

NXXl.NYYl.. NUMBER  OF  NODES  ALONG  X.Y-DIRECTIONS 

NYY . NUMBER  OF  DIVISIONS  BETWEEN  NODES  IN  Y-DIRECTION 

X( I). Yd).. COORDINATES  OF  THE  ITH  NODE 


Mesh  of  Quadrilateral  Elements  with  Four.  Eight,  or  Nine  nodes 


100  NEXl-NX+1 
NEYl-NY+1 
NXX-IEL*NX 
NYY-IEL*NY 
NXXl-NXX+1 


ooo  ooo  r>or> 


NVyi-NYY+l 

NEf^NX*NY 

NNM-NXX1*NYY1-(IEL-1)*NX*NY 

KO-0 

IF  (NPE.EQ.9)  THEN 
NNM-NXX1*NYY1 
KO-1 
END  IF 
C 

C  Generate  element  connectivity  array  *N0D(I,J>‘  of  first  element 
C 

N0D(1.1)-1 

N00(1.2)-1EL+1 

N00(1.3)-NXX1+(IEL-1)*NEX1+IEL+1 
IF  (NPE.EQ.9)  N0D(1.3)-4*NX+5 
N0D(1,4)-N0D(1.3)-IEL 
IF  (NPE.GT.4)  THEN 
N0D(1.5)-2 

N0D{1.6)-NXXl+(NPE-6) 

N0D(1.7)-N00(1.3)-1 
N00(1.8)-NXX1+1 
IF  (NPE.EQ.9)  NOO(1.9)-NXXl+2 
END  IF 


For  more  than  1  element  in  the  y-di recti  on 

200  IF  (NY.GT.l)  THEN 
M-1 

00  220  N-2.NY 
L-(N-1)*NX+1 
00  210  I-l.NPE 

210  N0D( L. I )-N00(M. I )+NXXl+( IEL-1 )*NEX1+K0*NX 

M-L 

220  CONTINUE 
END  IF 


For  more  than  1  element  in  the  x-direction 

230  IF  (NX.GT.l)  THEN 
00  260  NI-2.NX 
00  240  I-l.NPE 
Kl-IEL 

IF(I.EQ.6.0R.I.EQ.8)K1-1+K0 
240  N00(NI,I)-N0D(NI'l,I)-»-Kl 

M-NI 

00  260  NJ-2.NY 
L-(NJ-1)*NX+NI 
DO  250  J-l.NPE 

250  N0D( L. J )-N00(M. J )+NXXl+( IEL-1 )*NEX1+K0*NX 

M-L 

260  CONTINUE 
END  IF 

Generate  the  nodal  coorinates  arrays  'X(I)'  and  'Y(I)' 

270  YC-0.0 
C 

C  For  4  or  8-noded  elements 


c 


124 


IF  (NPE.EQ.9)  GOTO  310 
00  300  NI-l.NETl 

I-(NXX1+(1EL-1)*NEX1)*(NI-1)+1 

J-(NI-l)'iEL+l 

X{I)-0.0 

Y(1)-YC 

DO  280  NJ-l.NXX 
I-I+l 

X(I)-X(I-1)+0X(NJ) 

Y(I)-YC 
280  CONTINUE 
C 

C  For  8-nocled  elements 
C 

IF  (NI.GT.NY.OR.IEL.EQ.l)  GOTO  300 
J-J+1 

YC-YC+DY(J-1) 

I-I+l 

X(I)-0.0 

Y{I)-YC 

DO  290  II-l.NX 
K-2*II-1 
I-I+l 

X(I)-X(I-1)+DX{K)+0X(K+1) 
Y(I)-YC 
290  CONTINUE 

300  YC-YC+OY(J) 

RETURN 

C 

C  For  9-nocled  elements 
C 

310  DO  330  NI-1,NYY1 
I-NXX1*(NI-1) 

XC-0.0 

DO  320  NJ-l.NXXl 
I-I+l 
X(I)-XC 
YxI)-YC 
XC-XC+DX(NJ) 

320  CONTINUE 

YC-YC+OY{NI) 

330  CONTINUE 
C 

RETURN 

END 
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